Bias in error estimation when using cross-validation for model selection

Shrunken centroids

This classifier, originally proposed by Tibshirani et al [5], is an extension of the nearest centroid classifier. In the nearest centroid classifier the training set is used to calculate the centroids ${\overline{x}}_1$ and ${\overline{x}}_2$ (mean expressions of the genes) for the two classes. A new sample is compared to the two centroids and classified according to the class of the nearest centroid. In the shrunken centroids method, a parameter Δ is used to shrink the class centroids towards the overall centroid after standardizing by the within class standard deviation. The centroid belonging to class k is brought closer to the overall centroid $\overline{x}$ (mean of samples of all classes pooled together) by

where s is the vector of pooled within class standard deviation for all genes, s₀ is the median of the elements of s, d _k is given by

and (...)₊ denotes the positive part of the quantity in the parenthesis, i.e. equal to the quantity if it is greater than zero, and zero otherwise. Thus genes which are not very differentially expressed will contribute less to the classification than genes that are more discriminating. The parameter Δ can be varied to vary the number of genes used.

In the original paper, 10-fold CV is used to obtain the CV error estimate for a particular choice of Δ. There is no objective guideline given for selecting Δ based only on the training data. The authors vary Δ and use the value that minimizes the CV error estimate on the training set and the error on the testing data simultaneously. Thus the test set is used in the selection of the classifier parameters, which is problematic.

Support Vector Machines

Peng et al. present a method for feature selection using genetic algorithms (GA) and recursive feature elimination (RFE) in combination with a support vector machine (SVM) classifier [6]. SVMs were introduced by Vapnik [8] as linear hyperplanes that separate data belonging to different classes while maximizing the margin, or the distance of the training samples to the linear separating hyperplane.

Denote the two classes by 1 and -1. For a sample x consisting of p measurements (e.g. gene expressions) the linear hyperplane classifier c(x) predicts the class according to

for a weight vector w = ⌊w₀ w₁ ... w _p ⌋ and the augmented sample vector $\widehat{x}$ obtained by appending sample x with a constant 1, i.e. $\widehat{x}$ = ⌊1, x₁, ..., x _p ⌋.

The margin of a sample x of class y ∈ {-1,1} is defined as y$\widehat{x}$w' and it plays a very important role in SVM. The margin of correctly classified samples is positive and that for misclassified samples is negative. The margin of the classifier is defined as the smallest margin of all the training samples. The SVM tries to find a w such that the margin is maximized while the norm of the weight vector, ⟨w, w⟩ is minimized. This is equivalent to minimizing the cost function

where

ξ _i (w) = (1 - y _i ${\widehat{x}}_i$w ^t )₊     (5)

and (...)₊ denotes the positive part of the quantity in the parenthesis, as above, i.e. equal to the quantity if it is greater than zero, and zero otherwise.

Since the capacity of the classifier increases with increasing norm of the weight vector, the parameter C also controls the tradeoff between the size of the margin and the capacity of the classifier.

Since the samples need to be represented only in the form of scalar products, this formulation can be extended to non-linear classifiers by the introduction of kernels. Kernels are the functional representation of scalar products in transformed space. It can be shown that such a transformation leaves the optimization problem unchanged except that the inner product ⟨x₁, x₂⟩ is replaced by the kernel K(x₁, x₂). In the case of very high (or infinite) dimensional transformed space, the kernel is usually easier to compute than doing the transformation followed by scalar product.

For Gaussian kernel SVM (also called a radial basis function kernel), the kernel is given by

K(x₁, x₂) = exp(-γ ||x₁ - x₂||²)     (6)

The spread of the kernel function is given by γ, which can be varied to adapt the kernel to the data. The larger the value of γ, the more peaked the corresponding transformations of the feature vectors are, and the higher the capacity of the classifier.

Peng et al use the Leave-One-Out-CV (LOOCV) error to tune the kernel parameters. LOOCV is a CV scheme where one sample is left out during each iteration. The average classification error obtained is an almost unbiased estimate of the true error. In (6), the training data is used to select an appropriate kernel (from linear, Gaussian and polynomial) and set of parameters that minimize the LOOCV error estimate. However no independent test set is used and only the final LOOCV error estimate on the training set is reported.

Implementation

For each simulation we generated at least 1000 sets of 40 samples, of which 20 belonged to class 1 and the remaining 20 to class 2. Each sample was a vector of 6000 features (synthetic gene expressions). For some of the cases we used "null" data sets where no gene is differentially expressed between the two classes. For each gene, the population mean expressions in both classes were the same, namely zero. Thus none of the features are truly discriminatory between the two classes. We also used data generated from a "non-null" distribution for validating the nested CV approach. In this case, instead of using a "null distribution" (i.e. no difference in gene expression between the two classes), we simulated differential expression by fixing 10 genes (out of 6000) to have a population mean differential expression of 1 between the two classes. Samples from one class were drawn from a mean zero, unit variance, multivariate Normal distribution while samples for the second class had 10 genes with mean 1, unit variance Normal and the rest 5990 genes were of mean zero, unit variance Normal distribution.

Shrunken centroids

This was implemented in MATLAB^™ (Ver. 6.5, The Mathworks). For this case we used the "null" dataset with no difference between the two classes. For each 40 sample dataset, we computed the 10-fold CV error for the Shrunken Centroid classifier for values of Δ between 0.01 and 1. The value of Δ with the minimum 10-fold CV error CV(Δ) was determined.

Δ^* = arg min (CV(Δ)

If two values of Δ had the same error, the larger value was chosen (this corresponded to smaller number of genes with non-zero weights).

To determine if the CV error estimate CV(Δ^*) is a unbiased estimate of the true error for a classifier built using all the 40 samples with parameter value Δ^*, we created a Shrunken centroid classifier using all 40 samples and parameter value of Δ^*. We call this the optimized Shrunken Centroids classifier. This classifier was used to predict the class of 20000 samples created independently using the "null" data distribution. Since these test samples were not part of the training set, the mean error on them will give us the true error TE(Δ^*).

This process was repeated for each simulated 40-sample dataset providing empirical distributions of CV(Δ^*)and TE(Δ^*).

Support Vector Machine

The same type of analysis was performed for the SVM case. The "null" data distribution was used to create the synthetic data. For computational efficiency, we do not consider the complete algorithm used in (6). Instead of recursive feature elimination (RFE) for feature selection, we used the two sample t-statistics and selected the three features with the largest absolute t-statistic. In lieu of genetic algorithms (GA) as the search strategy, we used a simplified algorithm with a fixed Gaussian kernel. The classifier parameters tuned were C (trade off parameter in Eq. 4) and γ (kernel parameter). It must be noted that reducing the set of parameters over which one optimizes, as we have done here, may potentially reduce the amount of bias obtained.

The Leave-One-Out-CV (LOOCV) was used to compute the CV error estimate for each point on a grid of C and γ values (C was varied from 2^-5 to 2¹⁵ and γ from 2^-15 to 2³). The CV error estimate was computed by leaving out each of the 40 samples in turn, selecting the 3 best features using the t-statistic on the remaining 39 samples and creating an SVM classifier using the fixed C and γ values. This classifier was used to predict the class of the left out sample. The average error on the 40 samples is a CV error estimate CV(C, γ) for the values of C and γ used. This was repeated for all of the grid values. We used the LIBSVM package developed by Chang and Lin [9] to develop the classifier.

Similar to the analysis on Shrunken Centroids, we find the value of parameters that minimize the CV error estimate

(C^*, γ ^*) = arg min (CV(C, γ))

To compute the true error, an SVM classifier with parameters C^* and γ ^* was built using all 40 samples, and the top 3 features based on all 40 samples. We call this the optimized SVM classifier. This classifier was used to predict the classes of a large test set of 20000 samples with the same distribution as the training set (i.e. the "null" distribution). The mean error on the test set gives us the true error TE(C^*, γ ^*).

Nested CV with shrunken centroids and SVM

We evaluated the nested CV approach for the Shrunken Centroids classifier with Δ optimized using10-fold CV (the optimized Shrunken Centroids classifier). The "null" data distribution was used to generate the datasets. Finding the optimum value Δ^* of the classifier parameter was done exactly like described above, as was finding the true error TE(Δ^*). The only thing that changes is how a new estimate of the true error is computed. Instead of using the CV error estimate CV(Δ^*) for the optimal Δ, we used the nested CV error estimate. The nested CV error estimate is computed this way. One sample (out of the 40) was left out and the Δ for the Shrunken Centroids classifier was selected on the remaining 39 samples by minimizing the 10-fold CV. This is done exactly as above except that the wrapper algorithm is restricted to 39 samples, instead of 40. The wrapper algorithm determines an optimal value of Δ using the 39 samples and then creates a classifier with the same 39 samples and the optimal value selected based on these samples. This classifier is used to predict the class of the left out sample. This was done for each of the 40 samples, left out in turn. The average error over the samples is the nested CV error estimate CV _nest (Δ^*) for the optimized Shrunken Centroid classifier. Repeating this for several synthetic datasets gives us the empirical distribution of CV _nest (Δ^*) and TE(Δ^*).

The nested CV approach was also evaluated for the optimized SVM classifier. Here we used the "non null" data distribution to create the training samples (40 samples) and the test samples (20000 samples). Nested CV is used to find the nested CV error estimate CV _nest (C^*, γ ^*) of the SVM classifier optimized for parameters C and γ. The procedure is exactly the same as that described for optimized Shrunken Centroids above, except that the classifier training algorithm being evaluated is the optimized SVM classifier algorithm detailed earlier.

Fig 1 shows the empirical distributions of CV(Δ^*), the CV error estimate for the optimal Δ and TE(Δ^*), the true error for the optimal Δ for the optimized Shrunken Centroid classifier. Fig 2 shows the distributions for CV(C^*, γ ^*) and TE(C^*, γ ^*) for the optimized SVM classifier. In both plots, the CV error estimate is centered over means that are distinctly smaller than the mean true errors. Since both these simulations were done with "null" data, the true errors are centered on 50% while the CV error estimates have a lower mean. Numerical results are shown in Table 1. Even though the mean true error is 50% (i.e. equal to randomly choosing the classes), the CV error estimate on the training set averages 37.8% for the optimized Shrunken Centroid classifier and 41.7% for the optimized SVM classifier. On more than one-fifth of the random training datasets, the bias is more than 20% for the classifiers.

Fig 3 shows the distributions for the nested CV error estimate CV _nest (Δ^*) and the true error TE(Δ^*) for the optimized Shrunken Centroids. We obtain an almost unbiased estimate of the true error. The mean nested CV error estimate is a slight overestimate of the true error (54.2% compared to 50.0%), since the classifier used in each nested CV iteration is based on 39 samples, while the classifier used on the test set is trained on 40 samples.

Fig 4 shows the distribution of the nested CV error estimate CV _nest (C^*, γ ^*) and the true error TE(C^*, γ ^*) for the optimized SVM algorithm on the "non null" data sets. We obtain a mean nested CV error estimate of 35.3% compared to the true error of 32.0%. The higher estimate of training error compared to test error can again be attributed to the lower number of samples being used (39 vs. 40) to create the classifier in the nested CV iteration.