The C1C2: a framework for simultaneous model selection and assessment.

BACKGROUND: There has been recent concern regarding the inability of predictive modeling approaches to generalize to new data. Some of the problems can be attributed to improper methods for model selection and assessment. Here, we have addressed this issue by introducing a novel and general framework, the C1C2, for simultaneous model selection and assessment. The framework relies on a partitioning of the data in order to separate model choice from model assessment in terms of used data. Since the number of conceivable models in general is vast, it was also of interest to investigate the employment of two automatic search methods, a genetic algorithm and a brute-force method, for model choice. As a demonstration, the C1C2 was applied to simulated and real-world datasets. A penalized linear model was assumed to reasonably approximate the true relation between the dependent and independent variables, thus reducing the model choice problem to a matter of variable selection and choice of penalizing parameter. We also studied the impact of assuming prior knowledge about the number of relevant variables on model choice and generalization error estimates. The results obtained with the C1C2 were compared to those obtained by employing repeated K-fold cross-validation for choosing and assessing a model.
RESULTS: The C1C2 framework performed well at finding the true model in terms of choosing the correct variable subset and producing reasonable choices for the penalizing parameter, even in situations when the independent variables were highly correlated and when the number of observations was less than the number of variables. The C1C2 framework was also found to give accurate estimates of the generalization error. Prior information about the number of important independent variables improved the variable subset choice but reduced the accuracy of generalization error estimates. Using the genetic algorithm worsened the model choice but not the generalization error estimates, compared to using the brute-force method. The results obtained with repeated K-fold cross-validation were similar to those produced by the C1C2 in terms of model choice, however a lower accuracy of the generalization error estimates was observed.
CONCLUSION: The C1C2 framework was demonstrated to work well for finding the true model within a penalized linear model class and accurately assess its generalization error, even for datasets with many highly correlated independent variables, a low observation-to-variable ratio, and model assumption deviations. A complete separation of the model choice and the model assessment in terms of data used for each task improves the estimates of the generalization error.

Background

A common task in computational biology/bioinformatics and computational chemistry/chemometrics (hereafter abbreviated BBCC) is to model a dependent variable from a set of independent variables; this gives insight into the workings of the process being modeled and enables prediction of future observations. Typical examples include analyzing potential drug activity through proteochemometrics and quantitative structure-activity relationship (QSAR) modeling [1-3], discovering gene regulatory binding-site modules [4], and predicting clinical outcomes of cancer from gene expression data [5]. However, recent articles have indicated that predictive modeling approaches have not fully fulfilled expectations for solving real problems. This issue has for instance been discussed in the fields of QSAR [6] and microarray gene expression data modeling [7,8]. While some of the problems may be attributed to incorrect use and interpretation of the models, others can be ascribed to improper model selection and assessment. Our aim is here to address the latter issue by introducing the C1C2, a general framework for model choice and assessment.

Let = {, } be a dataset, where = (1, ..., )' is an n × pmatrix whose ith row, , is the value of a p-vector of independent variables associated with y, the ith row of the n × 1 matrix, = (y1, ..., y)'. A statistical model, can be used to characterize the relation between and . In general, given the dataset , must be chosen from a set of models, according to some criterion (typically the minimization of a loss function). The most common way to select from in BBCC is to use K-fold cross-validation; that is, the dataset is split into K mutually exclusive subsets, 1,..., ,..., , of approximately equal size and () is picked to minimize the function:

where L is a loss function and -denotes that the kth subset was excluded from during the model fitting; is the number of observations in subset D; and m = 1,..., M, where M is the number of models in . This model selection method may seem straightforward and intuitive, however it neglects the fact that all the data at hand is used to make the model choice. Thus, we no longer have an independent testset to assess the chosen model by. The result is that, typically, C0() underestimates the generalization error (see for instance [9]), defined as the expected prediction error over an independent test sample. This problem has been highlighted in relatively recent works [9,10], but was noted initially in 1974 [11]. To obtain a more accurate generalization error estimate, the model selection process must be separated from the model assessment in terms of the data that is used. Ideally, if data were abundant and easily produced, we would set aside a large test dataset and use it to assess – but not to choose! – the model , and subsequent model refinements could be assessed with new, unseen data. In practice, this is however often impossible since BBCC data is typically scarce, and expensive to produce. The luxury of large independent testsets can thus rarely be afforded. To tackle this problem, Freyhult et al. [9] suggested using a K-fold cross-validatory assessment of an H-fold cross-validatory model choice, , as a way of simultaneously choosing and assessing its performance; thereby separating the model selection from its assessment. The model is assessed by the function

where (-) is the cross-validatory choice of based on -; that is, the model that minimizes the function

where -denotes the dataset with the kth and hth subsets omitted. In the present work we build on and generalize this idea into the C1C2 framework. In general, the number of models in is huge, thus it is unfeasible to go through even a small subset of them manually. Hence, for a framework such as the C1C2 to be useful in practice, automated methods for searching the model space are necessary; in this sense the C1C2 is similar to the automatic modelling approaches taken in for instance [12-14]. Here, the specific use of the C1C2 is demonstrated by applying it together with two search methods to simulated and real-world datasets. The results are compared to those obtained by employing the function (1) for model selection and assessment. In the interest of clarity, we have restricted our attention to the study of model choice and assessment within a linear model class, (defined below) for a quadratic loss function. We discuss the results of the demonstrations, the pros and cons of the generality of the C1C2, and set out some directions for further research.

Results

Algorithm

Let C1, C2 ∈ = {C1,..., C}, where is a set of model assessment criteria and C1, C2 represent two specific criteria (i.e. C1 = C, C2 = C, i, j = 1,..., J). Further, let S ∈ , where is a set of search methods; let L ∈ , where is a set of loss functions; let denote a sequence of data processing steps (e.g. mean-centering, transformations, whitening, etc) and let ' contain the results of applied to -(the roles of and ' are exemplified in the discussion following the pseudococe below); let R be a positive integer and K an integer between 1 and n, where n is the number of observations. The C1C2 procedure is outlined with the following pseudocode:

Initiate , , L, C1, C2, R, and K.

for (r in 1,..., R) {

a. Partition data, = {}.

for (k in 1,..., K) {

b. Apply to -. Save results in '.

c. Search using the data -and C2 as objective function. Assume is found to maximize (or minimize) C2. Save .

d. Apply to using '.

e. Assess using C1 and . Save assessment result.

}

}

The data partitioning in (a) separates data for the model choice from data for the model assessment. Note that the partitioning is dependent on the choice of C1 and does not necessarily need to be done in a cross-validation fashion. For instance, the choice C1 = ".632 estimator" [15,16], partitions the data by independently sampling n rows from with replacements and lets the observations not included among the sampled observations constitute the test set. The output from the C1C2 is also dependent on the choice of C1; for example, the choice C1 = C2 = Bayesian Information Criterion (BIC, see [17] and Methods) would not give a direct estimation of the generalization error, but rather an assessment of model overfitting. To clarify the roles of and , we give the following example: Let only contain a processing step that scales to unit variance. In (b) is applied to -and the standard deviation of each column of -is saved in . In (d), is applied to , that is, the columns in are scaled using the standard deviations calculated in (b). This treatment of ensures that indeed constitutes an independent testset. The 'for loop' over r is introduced to enable calculation of confidence intervals for estimates and, by averaging the estimates over R repetitions, it permits reduction of the variance in parameter and error (or overfitting) estimates by a factor of 1/R.

Figure 1 gives a graphical gives a graphical view of the C1C2 framework. We emphasize that the generality of the framework allows C1, C2, and S to be chosen to fit the problem at hand. Adequate choices of C1, C2, and S make the model selection and assessment more accurate and faster, which we will discuss below.

Figure 1

The . The data partitioning in step (a) in the pseudocode separates the model choice from its assessment, which is highlighted in purple in the figure. The left side of the figure relates to steps (b) to (d) in the pseudocode, and the right side to step (e); i.e. the left side relates to choosing the model and saving the parameter estimates, and the right side to assessing the model and saving the assessment results.

Datasets used in the testing

Both the simulated and the real data used for evaluating a new method or algorithm should reflect typical dataset properties found in real-world application domains. Examples of such properties in BBCC are multicollinearity, a large number of independent variables relative to the number of observations, and binary and categorical independent variables.

Simulated data

We simulated datasets as follows:

Let Δ = (δ), where , represents a subset of , and let β(Δ) = (β(δ)), where , is a vector of regression coefficients. The data matrix, was sampled from a twenty-dimensional multivariate normal distribution. Thus, ~ N20(020, Σ20), i = 1,..., n, where 020 is a twenty-dimensional vector of zeroes and Σ20 is either the 20 identity matrix or a covariance matrix, 20, estimated from a real in-house QSAR dataset that originated from HIV protease inhibitors. The HIV QSAR dataset contains highly correlated independent variables resulting in an 20 with many large absolute values in the off-diagonal elements. Three δin Δ were chosen to be nonzero and equal to one; the positions were chosen at random to be 11, 14, and 18 (but remained fixed throughout the experiment for evaluation purposes). The corresponding regression coefficients, β11(δ11), β14(δ14), and β18(δ18), were, for simplicity, all set equal to 1. The variables 11, 14, and 18 were slightly correlated with an estimated covariance matrix .

Datasets were generated assuming that followed a linear model according to: = β(Δ) + ε, where ε~N(0,1.5). Four datasets were simulated in order to evaluate the C1C2s performance in settings where n p, and where the observations were sampled from an orthogonal multivariate normal distribution or not, according to the following schema:

1. n = 15, Σ20 = 20

2. n = 200, Σ20 = 20

3. n = 15, Σ20 = 20

4. n = 200, Σ20 = 20

The simulated datasets are available in CSV format from Additional files 1, 2, 3, 4.

The Selwood dataset

This is a real dataset, made available from a website [18] and originally published in 1990 [19]. It is a widely studied dataset (see [20,21] and references therein). It contains one dependent variable, 53 independent variables, and 31 observations. The 53 independent variables correspond to numerical descriptions of molecules (antifilarial antimycin analogues) designed to capture their physicochemical properties. The dependent variable is the in vitro antifilarial activity of the molecules. This dataset exhibits extremely strong correlations between the independent variables and contains real valued, binary and categorical independent variables. It is known from previous studies that this dataset contains nonlinearities, but that decent models can be found using linear methods.

Testing

To demonstrate the use of the C1C2, it was applied to the simulated and real-world datasets described above (hereafter referred to as "the datasets"). Below we describe the choices for R, K, , , L, C1, C2, and S and the motivation for each selection.

Choice of R and K

The larger the choice of R, the higher the accuracy in the estimate of the generalization error; the choice of R is thus constrained by time and is dependent on the size of the dataset and the computational complexity of the choices of , , L, C1, C2, and S. The choice of K is a trade-off between bias and variance; the larger the K, the more variance and the less bias in the estimates of the generalization error [22]. R was here set to 12 and K to 5.

Choice of

We make the assumption that a normal linear model forms a reasonable approximation of the data. This model is given by: y= β+ ε, where the subscripts denote the number of rows in a matrix, βare regression coefficients, and ε~N(0,σ2 ). Further, because n