Unbiased bootstrap error estimation for linear discriminant analysis

Convex bootstrap error estimation is a popular tool for classifier error estimation in gene expression studies. A basic question is how to determine the weight for the convex combination between the basic bootstrap estimator and the resubstitution estimator such that the resulting estimator is unbiased at finite sample sizes. The well-known 0.632 bootstrap error estimator uses asymptotic arguments to propose a fixed 0.632 weight, whereas the more recent 0.632+ bootstrap error estimator attempts to set the weight adaptively. In this paper, we study the finite sample problem in the case of linear discriminant analysis under Gaussian populations. We derive exact expressions for the weight that guarantee unbiasedness of the convex bootstrap error estimator in the univariate and multivariate cases, without making asymptotic simplifications. Using exact computation in the univariate case and an accurate approximation in the multivariate case, we obtain the required weight and show that it can deviate significantly from the constant 0.632 weight, depending on the sample size and Bayes error for the problem. The methodology is illustrated by application on data from a well-known cancer classification study. Electronic supplementary material The online version of this article (doi:10.1186/s13637-014-0015-0) contains supplementary material, which is available to authorized users.


Introduction
The bootstrap method [1][2][3][4][5][6][7] has been used in a wide range of statistical problems. The asymptotic behavior of bootstrap has been studied [8][9][10][11], while small-sample properties have been studied under simplifying assumptions, such as considering the estimator based on all possible bootstrap samples (the 'complete' bootstrap) [12][13][14]. The small-sample properties of the usual bootstrap are not well understood, in particular when it comes to estimating the error rates of classification rules [15,16].
There has been, on the other hand, interest in the application of bootstrap to error estimation in classification problems and, in particular, gene expression classification studies [17][18][19][20]. Of particular interest is the issue of classifier error estimation [21,22]. Bootstrap methods have generally been shown to outperform more traditional error estimation techniques, such as resubstitution and cross-validation, in terms of root-mean-square (RMS) error [4,5,7,[23][24][25][26][27][28][29][30][31][32][33][34][35]. Bootstrap error estimation is typically performed via a convex combination of the (generally) pessimistic basic bootstrap estimator, known as the zero bootstrap, and the (generally) optimistic resubstitution estimator. A basic problem is how to choose the weight that yields an unbiased estimator.
The problem of unbiased convex error estimation was previously considered in [36][37][38] for a convex combination of resubstitution and cross-validation estimators, and in [4,7,23] for a combination between resubstitution and the basic bootstrap estimator. In the former case, a fixed suboptimal weight of 0.5 was proposed in [36,38], while an asymptotic analysis to find the optimal weight was provided in [37]. In the latter case, our case of interest, a fixed suboptimal weight of 0.632 was proposed in [4], leading to the well-known 0.632 bootstrap estimator, while in [7], a suboptimal weight is computed by means of a samplebased procedure, which attempts to counterbalance the effect of overfitting on the bias, leading to the so-called 0.632+ bootstrap error estimator; the problem of finding the optimal weight for finite sample cases was addressed via a numerical approach in [23].
Here, we determine the optimal weight for finite sample cases analytically, in the case of linear discriminant analysis under Gaussian populations. In the univariate case, no http://bsb.eurasipjournals.com/content/2014/1/15 other assumptions are made. In the multivariate case, it is assumed that the populations are homoskedastic and that the common covariance matrix is known and used in the discriminant. In either case, no simplifications are introduced to the bootstrap error estimator; it is the usual one, based on a finite number of random bootstrap samples.
The analysis in this paper follows in the steps of previous papers that have provided analytical representations for the moments of error-estimator distributions [39,40]. In the univariate case, exact expressions are given for the expectation of the zero bootstrap error estimator, in the general heteroskedastic (general-variance) Gaussian case. By using similar expressions for the expected true and resubstitution error [39], this allows the exact calculation of the required weight. In the multivariate case, the expectation of the zero bootstrap error estimator is expressed as a probability involving the ratio of two noncentral chisquare variables, in the homoskedastic Gaussian case, assuming that the true common covariance matrix is used in the discriminant. The resulting expression is exact but necessitates approximation for its numerical computation. This is done in this paper via the Imhof-Pearson three-moment method, which is accurate in small-sample cases [41]. Use of similar expressions for the expected true and resubstitution error [40] then allows the exact calculation of the required weight.
In the homoskedastic case, the required weight for unbiasedness is shown to be a function only of the Bayes error and sample size. Accordingly, plots and tables of the required weight for varying values of Bayes error and sample size are presented; if the Bayes error can be estimated for a problem, this provides a way to obtain the optimal weight to use. In the univariate case, it was observed that as the sample size increases, the optimal weight settles on an asymptotic value of around 0.675, thus slightly over the heuristic value 0.632; by contrast, in the multivariate case (d = 2), the asymptotic value appears to be strongly dependent on the Bayes error, being as a rule significantly smaller than 0.632, except for very small Bayes error. This paper is organized as follows. The 'Bootstrap classification' section defines linear discriminant analysis as well as its application under bootstrap sampling. The 'Bootstrap error estimation' section reviews convex bootstrap error estimation. The 'Unbiased bootstrap error estimation' section contains the main theoretical results in the paper, providing the analytical expressions for the computation of the required convex bootstrap weight in the univariate and multivariate cases. The 'Gene expression classification example' section contains a demonstration of the usage of the optimal weight in bootstrap error estimation using data from the breast cancer classification study in [42,43]. Lastly, the 'Conclusions' section contains a summary and concluding remarks.
All the proofs are presented in the Appendix.

Bootstrap classification
Classification involves a predictor vector X ∈ R d , also known as a feature vector, which represents an individual from one of two populations 0 and 1 (we consider here only this binary classification problem). The classification problem is to assign X correctly to its population of origin. The populations are coded into a discrete label Y ∈ {0, 1}. Therefore, given a feature vector X, classification attempts to predict the corresponding value of the label Y. We assume that there is a joint feature-label distribution F XY for the pair (X, Y ) characterizing the classification problem. In particular, it determines the probabilities c 0 = P(X ∈ 0 ) = P(Y = 0) and c 1 = P(X ∈ 1 ) = P(Y = 1), which are called the prior probabilities. Given a fixed sample size n, the sample data is an The population-specific sample sizes are given by n 0 = n i=1 I Y i =0 and n 1 = n i=1 I Y i =1 = n − n 0 , which are random variables, with n 0 ∼ Binomial(n, c 0 ) and n 1 ∼ Binomial(n, c 1 ). When we need to emphasize that n 0 and n 1 are random variables, we will use capital letters N 0 and N 1 , respectively. This sampling design, which is the most commonly found one in contemporary pattern recognition, is known as mixture sampling [44].
A classification rule n is used to map the training data S n into a designed classifier ψ n = n (S n ), where ψ n is a function taking on values in the set {0, 1}, such that X is assigned to population 0 or 1 according to whether ψ n (X) = 0 or 1, respectively. The classification error rate ε n of classifier ψ n is the probability that the assignment is erroneous: where (X, Y ) is an independent test point and ε i n = P(ψ n (X) = 1 − i | Y = i) is the error rate specific to population i , for i = 0, 1. Since the training set S n is random, ε n is a random variable, with expected classification error rate E[ε n ]; this gives the average performance over all possible training sets S n , for fixed sample size n.
Linear discriminant analysis (LDA) employs Anderson's W discriminant [45], which is defined as follows: wherê are the sample means relative to each population, and is a matrix, which can be either (1) the true common covariance matrix of the populations, assuming it is known (this is the approach followed, for example, in [39,40,46]), or (2) the sample covariance matrix based on the pooled sample S n , which leads to the general LDA case. In this paper, we will assume case (1) throughout. The corresponding LDA classifier is given by that is, the sign of W (X) determines the classification of X.
A bootstrap sample S * n contains n instances drawn uniformly, with replacement, from S n . Hence, some of the instances in S n may appear multiple times in S * n , whereas others may not appear at all. Let C be a vector of size n, where the ith component C(i) equals the number of appearances in S * n of the ith instance in S n . The vector C will be referred to as a bootstrap vector.
For a given S n , the vector C uniquely determines a bootstrap sample S * n , which we denote by S C n . Note that the original sample itself is included: if C = (1, . . . , 1) def = 1 n , then S C n = S n , since each original instance appears once in the bootstrap sample. Note also that the number of distinct bootstrap samples, i.e., values for C, is equal to 2n−1 n ; even for small n, this is a large number. For example, the total number of possible bootstrap samples of size n = 20 is larger than 6.8 × 10 10 .
The vector C has a multinomial distribution with parameters (n, 1/n, . . . , 1/n), Starting from a classification rule n , one may design a classifier ψ C n = n (S C ) on a bootstrap training set S C . Its classification error ε C n is given as in (1) is the error rate specific to population i , for i = 0, 1. In this paper, we apply this scheme to the LDA classification rule defined previously. Notice the distinction between a bootstrap LDA classifier and a 'bagged' (bootstrapaggregated) LDA classifier [47,48]; these correspond to distinct classification rules. The bootstrap LDA classifier is employed here as an auxiliary tool to analyze the problem of unbiased bootstrap error estimation for the plain LDA classifier.

Bootstrap error estimation
Since the feature-label distribution is typically unknown, the classification error rate ε n has to be estimated by a sample-based statisticε n , commonly referred to as an error estimator. Data in practice are often limited, and the training sample S n has to be used for both designing the classifier ψ n and as the basis for the error estimatorε n . The simplest and fastest way to estimate the error of a designed classifier ψ n is to compute its error on the sample data itself: This resubstitution estimator, or apparent error, is often optimistically biased, that is, it is often the case that Bias ε r n = E ε r n −E[ε n ] < 0, though this is not always so. The bias tends to worsen with more complex classification rules [49].
The basic bootstrap error estimator is the zero bootstrap error estimator [4], which is introduced next. Given the training data S n , B bootstrap samples are randomly drawn from it. Denote the corresponding (random) bootstrap vectors by {C 1 , . . . , C B }. The zero bootstrap error estimator is defined as the average error committed by the B bootstrap classifiers on sample points that do not appear in the bootstrap samples: where n(C) is the number of zeros in C.
The bootstrap zero estimator tends to be pessimistically biased, since the amount of distinct training instances available for designing the classifier is on average (1 − e −1 )n ≈ 0.632n < n. Pessimistic bias in an error estimator can be mitigated by forming a convex combination with an optimistic error estimator [23]. In the case of bootstrap error estimation, the standard approach is to form a convex combination of the zero bootstrap with resubstitution, Selecting the appropriate weight w = w * leads to an unbiased error estimator, In [4], the weight w is heuristically set to w = 0.632 to reflect the average ratio of original training instances that appear in a bootstrap sample. This is known as the .632 bootstrap estimator which has been heavily employed in the machine learning field.

Unbiased bootstrap error estimation
The 0.632 bootstrap error estimator reviewed in the previous section is not guaranteed to be unbiased. In this http://bsb.eurasipjournals.com/content/2014/1/15 section, we will examine the necessary conditions for setting the weight w = w * in (8) to achieve unbiasedness. We will then particularize the analysis to the Gaussian linear discriminant case, where exact expressions for w * will be derived, both in the univariate and multivariate cases. The bias of the convex estimator in (8) is given by Setting this to zero yields the exact weight that produces an unbiased error estimator. Now, applying expectation on both sides of (7) produces where p(C) is given by (5) and the sum is taken over all possible values of C (an efficient procedure for listing all multinomial vectors is provided by the NEXCOM routine given in [50], Chapter 5). Equations (11) and (12) allow the computation of the weight w * given the knowledge of E[ ε n ], E ε r n , and E ε C n | C . We will present next exact formulas for these expectations in the case of the LDA classification rule under Gaussian populations.

Univariate case
In the univariate case, the common variance term cancels and the W statistic and LDA classifier become greatly simplified, with The following functions will be useful.
where Z is a zero-mean, unit-variance Gaussian random variable, and Z 1 , Z 2 are zero-mean, unit-variance random variables that are jointly Gaussian distributed, with correlation coefficient ρ.
Assume that population i is distributed as Under these conditions, John obtained in [39] an exact expression for the expectation of the true classification error for fixed sample sizes n 0 and n 1 (this is known as separate sampling [44]). John's result can be written as follows: where The corresponding result for E[ ε 1 n | N 0 = n 0 ] is obtained by simply interchanging all indices 0 and 1 in the previous expressions. The expected error rate can then be found by using conditioning and Equation (1): where P(N 0 = n 0 ) = n n 0 c n 0 0 c n 1 1 .
As for resubstitution, Hills provided in [51] exact expressions for the expected error for fixed n 0 and n 1 . However, his expression applies only to the case σ 0 = σ 1 . Theorem 3 in [52] provides a generalization of this result to the case of populations of unequal variances. First, note thatε r n = n 0 nε wherê are the apparent error rates specific to class 0 and 1, respectively. The result in [52] can be written as The corresponding result for E[ε r,1 n | N 0 = n 0 ] is obtained by interchanging all indices 0 and 1. The expected resubstitution error rate can then be found by using conditioning and Equation (18): Finally, let us consider the expected bootstrap error. Given C, the bootstrap LDA classifier is obtained by replacingμ i byμ C i , i = 0, 1, in (13): wherê are bootstrap sample means. Now, note that with N 0 = n 0 fixed, the training data labels Y i , i = 1, . . . , n, are no longer random. Since all classification rules of interest are invariant to reordering of the training data, we can, without loss of generality, reorder the sample points so that Y i = 0 for i = 1, . . . , n 0 , and Y 1 = 1 for i = n 0 + 1, . . . , n. Let the same reordering be applied to a given bootstrap vector C. The next theorem extends John's result to the classification error of the bootstrapped LDA classification rule defined by (23 where e = μ 1 − μ 0 The corresponding result for E ε C,1 n | N 0 = n 0 , C is obtained by interchanging all indices 0 and 1. Proof. See the Appendix. It is easy to check that the result in Theorem 1 reduces to the one in (14) and (15) when C = 1 n . Following (16), we can then write The expected bootstrap error rate E[ε boot n ] can now be computed via (12).
In the special case σ 0 = σ 1 = σ (homoskedasticity), it follows easily from the previous expressions that E[ ε n ], E[ε r n ], and E[ε boot n ] depend only on the sample size n and on the Mahalanobis distance between the populations δ = |μ 1 − μ 0 |/σ , and therefore so does the weight w * , through (11). Since the optimal (Bayes) classification error in this case is ε * = (−δ/2), there is a one-to-one correspondence between Bayes error and the Mahalanobis distance. Therefore, in the homoskedastic case, the weight w * is a function only of the Bayes error ε * and the sample size n. http://bsb.eurasipjournals.com/content/2014/1/15 We find that this value of M is large enough to obtain an accurate approximation. All other quantities are computed exactly, as described previously. One can see in Figure 1a that w * varies wildly and can be very far from the heuristic 0.632 weight; however, as the sample size increases, w * appears to settle around an asymptotic fixed value. This asymptotic value is approximately 0.675, being thus slightly larger than 0.632. In addition, Figure 1b allows one to see that convergence to the asymptotic value is faster for smaller Bayes errors. These facts help explain the good performance of the original convex 0.632 bootstrap error http://bsb.eurasipjournals.com/content/2014/1/15

Multivariate case
Assume that population i is distributed as a multivariate Gaussian N(μ i , ), for i = 0, 1. Under these conditions, John obtained in [39] an exact expression for the expectation of the error of the LDA classification rule, defined by (2) to (4), for the case where N 0 = n 0 is fixed. This result is stated by Moran in [40] as follows: where W 1 and W 2 are independently distributed as noncentral chi-square variables with d degrees of freedom http://bsb.eurasipjournals.com/content/2014/1/15 (d being the dimensionality) and noncentrality parameters λ 1 and λ 2 , with where is the squared Mahalanobis distance between the populations. The corresponding result for E[ ε 1 n | N 0 = n 0 ] is obtained by interchanging n 0 and n 1 . The expected true error rate can then be found by using (16).
Moran also provided the following expression for the expectation of the resubstitution error estimator in the multivariate case, for fixed N 0 = n 0 [40]: where W 3 and W 4 are independently distributed as noncentral chi-square variables with d degrees of freedom and noncentrality parameters λ 3 and λ 4 , with The corresponding result for E[ε r,1 n ] is obtained by interchanging n 0 and n 1 . The expected resubstitution error rate can then be found by using (22).
The bootstrap LDA classifier in the multivariate case is given by whereμ C 0 andμ C 1 are defined in (24). The next theorem generalizes John's result for the multivariate classification error to the case of the bootstrapped LDA classification rule. N(μ i , ), for i = 0, 1. Then, the expected error rate of the bootstrap LDA classification rule defined by (33) is given by

Theorem 2. Assume that population i is distributed as
where W 5 and W 6 are independently distributed as noncentral chi-square variables with d degrees of freedom and noncentrality parameters λ 5 and λ 6 , with where s 0 and s 1 are defined in (27). The corresponding result for E[ ε C,1 n | N 0 = n 0 , C] is obtained by interchanging s 0 and s 1 .
Proof. See the Appendix. It is easy to check that the result in Theorem 2 reduces to the one in (29) and (30) when C = 1 n .
An issue that arises in the multivariate case is the computation of the probabilities in (29), (31), and (34). This computation is very difficult since it involves the ratio of noncentral chi-square random variables, which has a doubly noncentral F distribution. Computation of this distribution is a hard problem. Moran proposes in [40] a complex procedure, based on work by Price [53], to compute this probability, which only applies to even dimensionality d. We employ a simpler procedure, namely, the Imhof-Pearson three-moment method, which is applicable to even and odd dimensionality [41]. This consists of approximating a noncentral χ 2 d (λ) random variable with a central χ 2 h random variable, by equating the first three moments of their distributions. This approach was also http://bsb.eurasipjournals.com/content/2014/1/15 employed in [52], where it was found to be very accurate. To fix ideas, we consider (29). The Imhof-Pearson three-moment approximation is given by where χ 2 h is a central chi-square random variable with h degrees of freedom, with and The approximation is valid only for c 3 > 0 [41]. If c 3 < 0, one uses the approximation where h and y are as in (37), and The same approximation method applies to (31) and (34) by substituting the appropriate values. As in the univariate case, the assumption of a common covariance matrix makes the expectations E[ ε n ], E[ε r n ], and E[ε boot n ] and thus also the weight w * , functions only of n and δ. Since ε * = (−δ/2), this means that the weight w * is a function only of the Bayes error ε * and the sample size n. Figure 2 and Table 2 display the value of w * computed with the previous expressions in this section, for several sample sizes and Bayes errors. As in the univariate case, E[ε boot n ] in (12) is approximated by a Monte Carlo procedure, with the same number M = 100 × n 2 of MC vectors. All other quantities are computed exactly, as described previously, save for the Imhof-Pearson approximation. We can see in Figure 2 that there is considerable variation in the value of w * and it can be far from the heuristic 0.632 weight; however, as the sample size increases, w * appears to settle around an asymptotic fixed value. In contrast to the univariate case, these asymptotic values here appear to be strongly dependent on the Bayes error and are significantly smaller than the heuristic 0.632 except for very small Bayes errors. As in the univariate case, convergence to the apparent asymptotic value is faster for smaller Bayes errors. These facts again help explain the good performance of the original convex 0.632 bootstrap error estimator for moderate sample sizes and small Bayes errors.

Gene expression classification example
Here we demonstrate the application of the previous theory in comparing the performance of the bootstrap error estimator using the optimal weight versus the use of the fixed w = 0.632 weight, using gene expression data from the well-known breast cancer classification study in [42], which analyzed expression profiles from 295 tumor specimens, divided into N 0 = 115 specimens belonging to the 'good-prognosis' population (class 1 here) and N 1 = 180 specimens belonging to the 'poor-prognosis' population (class 0).
Our experiment was set up in the following way. We selected two genes among the previously published 70gene prognosis profile [43]. These genes were selected for their approximate homoskedastic Gaussian distributions (see Figure 3). Since the real prior probabilities c 0 and c 1 for the good-and poor-prognosis populations are unknown, we assumed three different scenarios corresponding to c 0 = 1/3, c 0 = 1/2, and c 0 = 2/3 and downsampled randomly one or the other set of specimens to obtain new sample sizes (90, 180), (115, 115), and (115, 68), respectively, so as to reflect the assumed prior probabilities. In each of the three cases, we then drew 2,000 random samples of size n = 30 from the pooled data, computed for each the true error, resubstitution, basic bootstrap, and convex bootstrap error rates. Bias and root-mean-square (RMS) error for each estimator were estimated by averaging over the 2,000 repetitions. We considered both the fixed 0.632 weight and the optimal weight prescribed by our analysis. For the latter, we estimated for each value of c 0 the Bayes error using the full data set and read off Table 2 the optimal weight corresponding to the estimated Bayes error and sample size n = 30. The results are displayed in Table 3. Despite the approximate nature of the results, given that the simulated training samples are not independent from each other, we can see that the bias and RMS were always smaller for the estimator using the optimal weight than using the fixed 0.632 weight (all bootstrap estimators vastly outperforming resubstitution). given the sample size and problem difficulty, but also offer insight into the choice of the 0.632 weight for the classic 0.632 bootstrap error estimator. It was observed that the required weight for unbiasedness can deviate significantly from the 0.632 weight, particularly in the multivariate case, where the required weight for unbiasedness appears to settle on an asymptotic value that is strongly dependent on the Bayes error, being as a rule smaller than 0.632. The results were illustrated by application to gene expression data from a well-known breast cancer study.

Figure 3
Data used in the gene expression experiment. The plot shows the optimal (linear) classifier superimposed on the sample for the genes OXCT and WISP1, from the breast cancer study in [42]. We can see that both populations are approximately Gaussian with equal dispersion. Bad prognosis = red. Good prognosis = blue.
The result then follows after some algebraic manipulation. By symmetry, to obtain E[ ε 1 C | C], one needs only to interchange all indices 0 and 1.

Proof of Theorem 2
Following the same technique used in [32], we write E[ ε 0 C | C] = P(ψ C n (X) = 1 | X ∈ 0 , C) where U = (s 0 + s 1 ) − 1  where ρ c is defined as in (35) and I denotes the identity matrix of dimension d. It follows that are independent noncentral chi-squared random variables with d degrees of freedom and noncentrality parameters λ 5 and λ 6 defined in (35). The result then follows from (62). Following along the same lines, one can show that E[ ε 1 C | C] is obtained by interchanging s 0 and s 1 in the result for E[ ε 0 C | C] (the details are omitted for brevity).