Computing AIC for black-box models using Generalised Degrees of Freedom: a comparison with cross-validation

03/09/2016
by   Severin Hauenstein, et al.
University of Freiburg
0

Generalised Degrees of Freedom (GDF), as defined by Ye (1998 JASA 93:120-131), represent the sensitivity of model fits to perturbations of the data. As such they can be computed for any statistical model, making it possible, in principle, to derive the number of parameters in machine-learning approaches. Defined originally for normally distributed data only, we here investigate the potential of this approach for Bernoulli-data. GDF-values for models of simulated and real data are compared to model complexity-estimates from cross-validation. Similarly, we computed GDF-based AICc for randomForest, neural networks and boosted regression trees and demonstrated its similarity to cross-validation. GDF-estimates for binary data were unstable and inconsistently sensitive to the number of data points perturbed simultaneously, while at the same time being extremely computer-intensive in their calculation. Repeated 10-fold cross-validation was more robust, based on fewer assumptions and faster to compute. Our findings suggest that the GDF-approach does not readily transfer to Bernoulli data and a wider range of regression approaches.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

09/05/2018

Cross validation residuals for generalised least squares and other correlated data models

Cross validation residuals are well known for the ordinary least squares...
11/15/2017

Accelerating Cross-Validation in Multinomial Logistic Regression with ℓ_1-Regularization

We develop an approximate formula for evaluating a cross-validation esti...
06/16/2020

On parametric tests of relativity with false degrees of freedom

General relativity can be tested by comparing the binary-inspiral signal...
10/21/2019

Generalised learning of time-series: Ornstein-Uhlenbeck processes

In machine learning, statistics, econometrics and statistical physics, k...
01/18/2021

The Violating Assumptions Series: Simulated demonstrations to illustrate how assumptions can affect statistical estimates

When teaching and discussing statistical assumptions, our focus is often...
08/22/2019

Efficient Cross-Validation of Echo State Networks

Echo State Networks (ESNs) are known for their fast and precise one-shot...
06/01/2018

Return of the Infinitesimal Jackknife

The error or variability of machine learning algorithms is often assesse...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In many scientific fields, statistical models have become a frequently used tool to approach research questions. Choosing the most appropriate model(s), i.e. the model(s) best supported by the data, however, can be difficult. Especially in ecological, sociological and psychological research, where data are often sparse while systems are complex, the evidence one particular statistical model may not be exclusive. The existence of alternative models, which fit the data with comparable goodness, but yield considerably different predictions is ubiquitous (e.g. Draper 1995, Madigan et al. 1995, Raftery 1996).

In ecological and evolutionary research, statistical procedures are dominated by an “information-theoretical approach” (Burnham & Anderson 1998, 2002), which essentially means that model fit is assessed by Akaike’s Information Criterion (defined as : Akaike 1973). In this field, the AIC has become the paradigmatic standard for model selection of likelihood-based models as well as for determination of model weights in model averaging (e.g. Diniz-Filho et al. 2008, Mundry 2011, Hegyi & Garamszegi 2011). In recent years, and particularly in the field of species distribution analyses, non-parametric, likelihood-free approaches (‘machine learning’) have become more prevalent and in comparisons typically show better predictive ability (e.g. Recknagel 2001, Elith et al. 2006, Olden et al. 2008, Elith & Leathwick 2009). However, these approaches do not allow the computation of an AIC, because many of such ‘black-box’ methods are neither likelihood-based, nor can one readily account for model complexity, as the number of parameters does not reflect the effective degrees of freedom (a phenomenon Elder (2003)

calls the “paradox of ensembles”). Thus, at the moment ecologists are faced with the dilemma of either following the AIC-paradigm, which essentially limits their toolbox to GLMs and GAMs, or use machine-learning tools and closing the AIC-door. From a statistical point of view, dropping an AIC-based approach to model selection and averaging is no loss, as an alternative approximation of the Kullback-Leibler distance is possible through cross-validation.

In this study, we explore a potential avenue to unite AIC and machine learning, based on the concept of Generalised Degrees of Freedom (GDF). We explore the computation of GDFs for Gaussian and Bernoulli-distributed data as plug-in estimates of the number of parameters in a model. We also compare such a GDF-AIC-based approach with cross-validation.

The paper is organised as follows. We first review the Generalised Degrees of Freedom concept and relate it to the degrees of freedom in linear models. Next, we briefly illustrate the relation between cross-validation and Kullback-Leibler divergence, as KL also underlies the derivation of the AIC. Through simulation and real data we then explore the computation of GDFs for a variety of modelling algorithms and its stability. The paper closes with a comparison of GDF-based AIC and cross-validation-derived deviance, and comments on computational efficiency of the different approaches and consequences for model averaging.

1.1 Generalised Degrees of Freedom

Generalised Degrees of Freedom (GDF), originally proposed by Ye (1998) and illustrated for machine learning by Elder (2003), can be used as a measure of model complexity through which we can make information theoretical approaches applicable to blackbox algorithms.

In order to understand the conceptual properties of this method, we can make use of a somewhat simpler version of degrees of freedom (df). For a linear model m, , df are computed as the trace of the hat (or projection) matrix H, with elements (e.g. Hastie et al. 2009, p. 153):

(1)

For a linear model this is the number of its independent parameters, i.e. the rank of the model. To expand this concept to other, non-parametric methods, making it in fact independent of the actual fitted model, the definition above provides the basis for a generalised version of degrees of freedom. Thus, and according to Ye (1998),

(2)

where is the fitted value. That is to say, a model is considered the more complex the more adaptively it fits the data, which is reflected in higher GDF. Or in the words of Ye (1998, p. 120): GDF quantifies the “sensitivity of each fitted value to perturbations in the corresponding observed value”. For additive error models (i.e. , with ), we can use the specific definition that

(3)

(Hastie et al. 2009, p. 233). Beyond the additive-error model more generally, however, we have to approximate eqn 2 in other ways. By fitting the model with perturbed and assessing the response in we can evaluate so that,

(4)

where is perturbed in such a way that, for normally distributed data, with

being small relatively to the variance of

.

In order to adapt GDF to binary data we need to reconsider this procedure. Since a perturbation cannot be achieved by adding small random noise to , the only possibility is to replace 0s by 1s and vice versa. A perturb-all-data-at-once approach (Elder 2003) as for Gaussian data is thus not feasible. We explore ways to perturb Bernoulli-distributed below.

The equivalence of GDF and in the linear model encourages us to use GDF as a plug-in estimator of the number of parameters in the AIC computation.

1.2 Cross-validation and a measure of model complexity

Cross-validation is a standard approach to quantify predictive performance of a model, automatically accounting for model complexity (by yielding lower fits for too simple and too complex models, e.g. Wood 2006, Hastie et al. 2009). Because each model is fitted repeatedly, cross-validation is computationally more expensive that the AIC, but the same problem arises for GDF: it requires many simulation to estimate it stably (see below).

We decided to use the log-likelihood as measure of fit for the cross-validation (Horne & Garton 2006), with the following reasoning. Let denote the model density of our data, y, and let denote the true density. Thus the Kullback-Leibler (KL) divergence (Kullback & Leibler 1951) is

(5)

AIC is an estimate of this values (in the process taking an expectation over the distribution of ). Now is independent of the model, and the same for all models under consideration (and hence gets dropped from the AIC). So the important part of the KL-divergence is

(6)

that is, the expected model log-likelihood, where the expectation is over the distribution of data not used in the estimation of . Expression 6 can be estimated by cross-validation:

(7)

where is the sum over K folds of the log likelihood of the test subset , given the trained model , with parameter estimates . To put this on the same scale as AIC we multiply by to obtain the cross-validated deviance.

With the assumption of AIC and (leave-one-out) cross-validation being asymptotically equivalent (Stone 1977) and given the definition of AIC, we argue that it should be possible to extract a measurement from that quantifies model complexity. Hence,

(8)

with representing (estimated) model complexity, the maximum log-likelihood of the original (non-cross-validated) model and N the number of data points. Embracing the small sample-size bias adjustment of AIC (Sugiura 1978, Hurvich & Tsai 1989) we get:

(9)

Thus, we can compute both a cross-validation-based deviance that should be equivalent to the AIC, , as well as a cross-validation alternative to GDF, based on the degree of overfitting of the original model (eqns. 8 and 9).

This approach to computing model complexity is not completely unlike that of the DIC (Spiegelhalter et al. 2002), where represents the effective number of parameters in the model and is computed as the mean of the deviances in an MCMC-analysis minus the deviance of the mean estimates: (Wood 2015). In eqn. 8 the likelihood estimate plays the role of , while the cross-validation estimates .

2 Implementing and evaluating the GDF-approach for normally and Bernoulli-distributed data

We analysed simulated and real data sets using five different statistical models. Then we computed their GDF, disturbing data points at a time, with different intensity (only for normal data), and for different values of .

2.1 Data: simulated and real

First, we evaluated the GDF-approach on normally distributed data following Elder (2003), deliberately using a relatively small data set: . The response was simulated as (, with -values of 10 and 10, respectively), with . (This simulation was repeated to control for possible influences of the data itself on the resulting GDF, but results were near-identical.)

We simulated binary data with (effective sample size ESS) and , (with -values of , 5, , 10, respectively), with .

Our real-world data comprised a fairly small data set (, ESS), showing the occurrence of sperm whales (Physeter macrocephalus) around Antarctica (collected during cetacean IDCR-DESS SOWER surveys), and a larger global occurrence data set (, ESS) for the red fox (Vulpes vulpes, provided by the International Union for Conservation of Nature (IUCN)). We pre-selected predictors as described in (Dormann & Kaschner 2010, Dormann et al. 2010), yielding the six and three (respectively) most important covariates. Since we are not interested in developing the best model, or comparing model performance, we did not strive to optimise the model settings for each data set.

2.2 Implementing GDF for normally distributed data

A robust computation of GDF is a little more problematic than eq. 4 suggests. Ye (1998)

proposed to fit a linear regression to repeated perturbations of each data point, i.e. to

over for each repeatedly perturbed data point , calculating GDF as the sum of the slopes across all data points. As is constant for all models, and is constant for the non-stochastic algorithms (GLM, GAM), the linear regression simplifies to over . For internally stochastic models we compute a mean as plugin for , and therefore apply the procedure also to randomForest, ANN and BRT.

Elder (2003) presents this so-called “horizontal” method (the perturbed and fitted are stored as different columns in two separate matrices, which are then regressed row-wise) as more robust compared to the “vertical” method (where for each perturbation a GDF is computed for each column in these matrices and later averaged across replicate perturbations). Convergence is poor when using the “vertical” method, so we restricted the computation of GDF to the “horizontal” method.

2.3 GDF for binary data

While normal data can be perturbed by adding normal noise (see next section), binary data cannot. The only option is to invert their value, i.e. or . Clearly, we cannot perturb many data points simultaneously this way, as that will dilute any actual signal in the data. However, for large data sets it is computationally expensive to perturb all data points individually, and repeatedly (to yield values for the horizontal method, see above). We thus varied the number of data points to invert, , to evaluate whether we can raise without biasing the GDF estimate.

2.4 How many data points to perturb simultaneously?

With the number of points to perturb, , the number of replicates for each GDF-computation () and the amplitude of disturbance (for the normal data), we have three tuning parameters when computing GDF.

First, we calculated GDF for the simulated data and the sperm whale dataset with taking up increasing values (from 1 to close to , or to ESS for binary data). Thus, random subsets of size

of the response variable

were perturbed, yielding . After each particular perturbation the model was re-fitted to . To gain insight about the variance of the computed GDF values we repeated this calculation 100 times (due to very high computational effort, the number of reruns had to be limited to 10 for both randomForest and BRT).

The perturbation for the normally distributed were also drawn from a normal distribution . We evaluated the sensitivity to this parameter by setting it to 0.125 and 0.5.

2.5 Modelling approaches

We analysed the data using Generalised Linear Model (GLM), Generalised Additive Model (GAM), randomForest, feed-forward Artificial Neural Network (ANN) and Boosted Regression Trees (BRT). For the GLM, GDF should be identical to the trace of the Hessian (and the rank of the model), hence the GLM serves as benchmark. For GAMs, different ways to compute the degrees of freedom of the model have been proposed (Wood 2006), while for the other three methods the intrinsic stochasticity of the algorithm (or in the case of ANN of the initial weights) may yield models with different GDFs each time.

All models were fitted using R (Team 2014) and packages gbm (for BRTs: Ridgeway et al. 2013, interaction depth=3, shrinkage=0.001, number of trees=3000, cv.folds=5), mgcv (for GAMs: Wood 2006, thin plate regression splines), nnet (for ANNs: Venables & Ripley 2002, size=7, decay=0.03 for simulated and decay=0.1 for real data, linout=TRUE for normal data), randomForest (Liaw & Wiener 2002) and core package stats (for GLMs, using also quadratic terms and first-order interactions). R-code for all simulation as well as data are available on https://github.com/biometry/GDF.

2.6 Computation of AIC and AIC-weights, from GDF and cross-validation

In addition to the number of parameters, the AIC-formula requires the likelihood of the data, given the model. As machine-learning algorithms may minimise a score function different from the likelihood, the result probably differs from a maximum likelihood estimate. To calculate the AIC, however, we have to assume that the distance minimised by non-likelihood methods is proportional to the likelihood, otherwise the AIC would not be a valid approximation of the KL-distance. No such assumption has to be made for cross-validation, as

here only serves as a measure of model performance. For the normal data, we compute the standard deviation of the model’s residuals as plug-in estimate of

.

For the binary data, we use the model fits as probability in the Bernoulli distribution to compute the likelihood of that model. We then calculated the AIC for all considered models based on their GDF-value. Due to the small sample sizes, we used AICc (Sugiura 1978, Hurvich & Tsai 1989):

(10)

We used (10-fold) cross-validation, maintaining prevalence in the case of binary data, to compute the log-likelihood of the cross-validation, yielding . To directly compare it to AICc, we multiplied with . The cross-validation automatically penalises for overfitting by making poorer predictions.

For model averaging, we computed model weights for each model , once for the GDF-based AICc and for the cross-validation log-likelihood, using the equation for Akaike-weights (Turkheimer et al. 2003, Burnham & Anderson 2002, p. 75):

(11)

where , taking the smallest AICc, i.e. the AICc of the best of the candidate models as ; n is the number of models to be averaged over.

The same idea can be applied to cross-validated log-likelihoods, so that

(12)

where , with being the largest cross-validated log-likelihood in the model set.

3 Results

3.1 GDF configuration analysis

For normally distributed data, increasing the number of points perturbed simultaneously typically slightly increased the variance of the Generalised Degrees of Freedom calculated for the model (Fig. 1 left column, GLM, GAM, but not for randomForest and ANN). For GLM and GAM, GDF computations yielded exactly the same value as the model’s rank (indicated by the dashed horizontal line).

For simulated Bernoulli data (Fig. 1 central column), we also observed an effect of the number of points perturbed on the actual GDF value, which decreased for GLM, GAM and ANN, but increased for BRTs with the number of data points perturbed. Several data points needed to be perturbed () to yield an approximately correct estimate. More worryingly, GDF depended non-linearly on the number of data points perturbed, with values varying by a factor of 2 for GAM and ANN. For GAM the sensitivity occurs because the smoothing parameter selection is sensitive to the quite severe information loss as more and more data are perturbed. GLM and BRT yielded more consistent but still systematically varying GDF-estimates.

The same pattern was observed for the sperm whale data (Fig. 1, right column). For the GLM, there was still some bias observable, also in the sperm whale-case study’s GLM. We attribute it to the fact that by perturbing the binary data we also alter the prevalence.

GaussianBernoulliSperm WhaleGLMGAMrandomForestANNBRTdata points perturbed simultaneously (k)

Figure 1: Models’ GDFs as a function of the numbers of data points perturbed (k) (ranked abscissa) for the simulated Gaussian (left), Bernoulli (centre) and sperm whale data (right column) and the five model types (rows). The filled red and black dot represent the mean GDF for the two replications at each level of k (in light-grey open dots and black open dots ). The dashed line in the first row       represents the number of model parameters of the GLM.

The two replicate simulations yielded consistent estimates, except in the case of the normal GAM and normal BRT. This suggests that both methods “fitted into the noise” specific to the data set, while the other methods did not. We did not observe this phenomenon with the binary data.

For randomForest, GDF-estimates centre around , meaning that perturbations of the data did not affect the model predictions. This is to some extent explicable by the stochastic nature of randomForest, giving different predictions when fitted to the same data. We interpret the value of as the perturbations creating less variability than the intrinsic stochasticity of this approach.

This is not a general feature of stochastic approaches, as in neural networks and boosted regression trees the intrinsic stochasticity seems to be much less influential, and both approaches yielded relatively consistent GDF-values. The actual GDF estimate of course depends on the settings of the methods and should not be interpreted as representative.

To compute the GDF for normally distribute data, we have to choose the strength of perturbation. The GDF-value is robust to this choice, unless many data points are disturbed (Fig. 2). Only for BRT does increasing the strength of perturbation lead to a consistent, but small, decrease in GDF-estimates, suggesting again that BRTs fit into the noise.

Estimates of GDF for randomForest and ANN (but not BRT) respond with decreasing variance to increasing the strength of perturbation. Increasing the intensity of perturbation seems to overrule their internal stochasticity.

GLMGAMrandomForestANNBRTperturbed data points at once (k)

Figure 2: Perturbation magnitude: Estimated GDF along an increasing numbers of perturbed data points for the simulated Gaussian data. For each there are 100 GDF replications computed with (), () and (), hence lower to higher perturbation magnitude. The dotted green line         represents the model’s rank (GLM) and the sum of estimated degrees of freedom (GAM). Note that x-axis is rank-transformed.

It seems clear from the results presented so far that no single best perturbation strength and optimal proportion of data to perturb exists for all methods. For the following analyses, we use, for normal data, perturbation and (except for GAM, where ); and for Bernoulli data (except for BRT and ANN, where ).

Figure 3: Development of mean Generalised Degrees of Freedom (GDF, left) and cross-validation log-likelihood-differences (, right) over 1000 replicates (250000 and 10000 model runs, respectively). GDF computation with data points perturbed per model run and 50 internal replicates for the Gaussian simulation data (), i.e. 250 model evaluations per replicate. derives from 10-fold cross-validation, i.e. 10 model evaluations per replicate. The dashed lines display the mean GDF and , respectively, of all replicates

one standard error.

3.2 Efficiency of GDF and cross-validation computations

Both GDF and cross-validation require multiple analyses. For the GDF, we need to run several perturbations, and possibly replicate this many times. For cross-validation, we may also want to repeatedly perform the 10-fold cross-validation itself to yield stable estimates. The mean GDF and CV-log-likelihood () over 1-1000 replicates are depicted in Fig. 3. With 100 runs, both estimates have stabilised, but 100 runs for GDF represent 25000 model evaluations (due to the perturbations and 50 internal replicates for a data set of ), while for 10-fold cross-validation these represent only 1000 model evaluations, making it 25 times less costly.

3.3 GDF-based AIC vs Cross-Validation

We analysed four data sets with the above settings for GDF and cross-validation, two simulated (Gaussian and Bernoulli) and two real-world data sets (sperm whale and red fox geographic distributions).

Both Generalised Degrees of Freedom (GDF) and cross-validation log-likelihood-differences () measure model complexity in an asymptotically equivalent way (section 1.2). For finite data sets both approaches yield identical rankings of model complexity but are rather different in absolute value (eqn. 1). Particularly the red fox-case study yields very low GDF-estimates for GLM, GAM and randomForest, while their cross-validation model complexity is much higher.

Model GDF GDF
Gaussian simulation Bernoulli simulation
GLM
GAM
randomForest
ANN
BRT
Sperm whale Red fox
GLM
GAM
randomForest
ANN
BRT
Table 1: Measures of model complexity (mean standard deviation): Generalised Degrees of Freedom (GDF) and derived from cross-validation (see section 1.2) for Gaussian and Bernoulli simulation data and real-world sperm whale and red fox distribution data. The true ranks for the GLMs are 15, 15, 7 and 8, respectively.

Model complexity is only one term in the AIC-formula (eqn. 8) and for large data sets the log-likelihood will dominate. To put the differences between the two measures of model complexity into perspective, we computed AIC for all data sets and compared it to the equivalent cross-validation deviance ().

Figure 4: Cross-validated deviance versus GDF-based AICc. Error bars display one standard error (sometimes too small to be visible). The identity line () indicates the equivalence of both measures. Colours represent the four datasets: Gaussian simulated (top-left, ) and Bernoulli simulated (top-centre, ); sperm whale (bottom-left, ) and red fox (bottom-centre, ) data. In the right panel per datum GDF and , respectively, i.e. divided by the number of data points.

Across the entire range of data sets analysed both approaches yield extremely similar results (Fig. 4 right). Within each data set, however, the pattern is more idiosyncratic, revealing a high sensitivity to low sample size (, i.e. all data sets except the red fox).

3.4 Gdf, and model weights

For all four datasets one modelling approach always substantially outperforms the others, making model averaging an academic exercise. We compare model weights (according to eqns 11 and 12) purely to quantify the differences in AIC and cross-validation deviance for model averaging. Only for the Bernoulli simulation is the difference noticeable (Table 2). Here GLM and randomForest share the model weight when quantified based on cross-validation, while for the GDF-approach GLM takes all the weight.

Model
Gaussian simulation Bernoulli simulation
GLM
GAM
rF
ANN
BRT
Sperm whale Red fox
GLM
GAM
rF
ANN
BRT

Table 2: Model weights: Comparison of Akaike weights and cross-validation weights (see eq. 12) for Gaussian and Bernoulli simulation data and real-life sperm whale and red fox abundance data. The “surprisingly” good performance of GLMs in the sperm whale case study remains unexplained, but is fully reproducible.

4 Discussion

So far, Ye’s (1998) Generalised Degrees of Freedom-concept did not attract much attention in the statistical literature, even though it builds on established principles and applies to machine learning, where model complexity is unknown. Shen & Huang (2006) have explored the perturbation approach to GDFs in the context of adaptive model selection for linear models and later extended it to linear mixed effect models (Zhang et al. 2012). The also extended it to the exponential family (including the Bernoulli distribution) and even to classification and regression trees (Shen et al. 2004). Our study differs in that the models considered are more diverse and internally include weighted averaging, which clearly poses a challenge to the GDF-algorithm.

4.1 GDF for normally distributed data

For normally distributed data our explorations demonstrate a low sensitivity to the intensity of perturbation used to compute GDFs. Furthermore, across the five modelling approaches employed here, GDF estimates are stable and constant for different numbers of data points perturbed simultaneously.

GDF-estimates were consistent with the rank of GLM models and in line with the estimated degrees of freedom reported by the GAM. For neural networks and boosted regression trees, GDF-estimates appear plausible, but cannot be compared with any self-reported values. Compared to the cross-validation method, GDF values are typically, but not consistently, lower by 10-30% (Table 1).

For randomForest, GDF-estimates were essentially centred on zero. It seems strange to find that an algorithm that uses hundreds of classification and regression trees internally to actually have no (or even negative) degrees of freedom. We expected a low value due to the averaging of submodels (called the ‘paradox of ensembles’ by Elder 2003), but not such a complete insensitivity to perturbations. (Within this study, SH and CFD independently reprogrammed our GDF-function to make sure that this was not due to a programming error.) As eqn. 4 shows, the perturbation of individual data point are compared to the change in the model expectation for this data point, and then summed over all data points. To yield a GDF of 0, the change in expectation (numerator) must be much smaller than the perturbation itself (denominator). This is possible when expectations are variable due to the stochastic nature of the algorithm. It seems that randomForest is much more variable than the other stochastic approaches of boosting and neural networks.

4.2 Bernoulli GDF

Changing the value of Bernoulli-data from 0 into 1 (or vice versa) is a stronger perturbation than adding a small amount of normal noise to Gaussian data. As our exploration has shown, the GDF for such Bernoulli-data is indeed much less well-behaved than for the normal data. Not only is the estimated GDF dependent on the number of data points perturbed, also is this dependence different for each modelling approach we used. This makes GDF-computation impractical for Bernoulli data. As a consequence, we did not attempt to extend GDFs in this way to other distributions, as in our perception only a general, distribution- and model-independent algorithm is desirable.

4.3 GDF vs model complexity from cross-validation

Cross-validation is typically used to get a non-optimistic assessment of model fit (e.g. Hawkins 2004). As we have shown, it can also be used to compute a measure of model complexity similar (in principle) to GDF (eqn. 8 and Table 1). Both express model complexity as the effective number of parameters fitted in the model. GDF and cross-validation-based model complexity estimator are largely similar, but may also differ substantially (Table 4, red fox case study). Since the ‘correct’ value for this estimator is unknown, we cannot tell which approach actually works better. Given our inability to choose the optimal number of data points to perturb (except for GLM), we prefer , which does not make any such assumption.

4.4 Remaining problems

To make the GDF approach more generally applicable, a new approach has to be found. The original idea of Ye (1998) is appealing, but not readily transferable in the way we had hoped.

Another problem, even for Gaussian data where this approach seems to be performing fine, is the high computational burden. GDF-estimation requires tens of thousands of model evaluations, giving it very limited appeal, except for small data sets and fast modelling approaches. Cross-validation, as alternative, is at least an order of magnitude faster, but still requires around 1000 evaluations. If the aim is to compute model weights for model averaging, no precise estimation of model complexity is needed and even the results of a single 10-fold cross-validation based on eqn. 8 can be used. It was beyond the scope of this study to develop an efficient cross-validation-based approach to compute degrees of freedom, but we clearly see this as a more promising way forward than GDF.

4.5 Alternatives to AIC

The selection of the most appropriate statistical model is most commonly based on Kullback-Leibler (KL) discrepancy (Kullback & Leibler 1951): a measure representing the distance between the true and an approximating model. Thus, we assume that a model , for which the distance to the true model is minimal, is the KL-best model. Yet, since KL-discrepancy is not observable, even if a true model existed, many statisticians have attempted to find a metric approximation (e.g. Burnham & Anderson 2002, Burnham et al. 1994). Akaike (1973), who proposed this measure as the basis for model selection in the first place, developed the AIC to get around the discussed problem.

The point of the cross-validated log-likelihood is that we do away with the approximation that yields the degrees of freedom term in the AIC, instead estimating the model-dependent part of the KL divergence directly. This approach is disadvantageous if AIC can be computed from a single model fit. But if the EDF terms for the AIC would require repeated model fits then there is no reason to use the AIC-approximation to the KL-divergence, rather than a more direct estimator. If leave-one-out cross-validation is too expensive, then we can leave out several, at the cost of some Monte-Carlo variability (resulting from the fact that averaging over all possible left out sets is generally impossible).

5 Conclusion

We have shown that the idea of using GDFs to extend information-theoretical measures of model fit (such as AIC) to non-likelihood models is burdened with large computational costs and yields variable results for different modelling approaches. Cross-validation is more variable than GDFs, but a more direct way to compute measures of model complexity, fit and weights (in a model averaging context). As cross-validation may, but need not, employ the likelihood fit to the hold-out, it appears more plausible for models that do not make likelihood assumptions. Thus, we recommend repeated ( times) 10-fold cross-validation to estimate any of the statistics under consideration.

Acknowledgements

This work was partly performed on the computational resource bwUniCluster funded by the Ministry of Science, Research and Arts and the Universities of the State of Baden-Württemberg, Germany, within the framework program bwHPC.

References

  • Akaike (1973) Akaike, H. (1973). Information theory and an extension of the maximum likelihood principle. In Second International Symposium on Information Theory (pp. 267–281).: Akademinai Kiado.
  • Burnham & Anderson (1998) Burnham, K. P. & Anderson, D. R. (1998). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Berlin: Springer.
  • Burnham & Anderson (2002) Burnham, K. P. & Anderson, D. R. (2002). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Berlin: Springer. 2nd edition.
  • Burnham et al. (1994) Burnham, K. P., Anderson, D. R., & White, G. C. (1994). Evaluation of the Kullback-Leibler discrepancy for model selection in open population capture-recapture models. Biometrical Journal, 36(3), 299–315.
  • Diniz-Filho et al. (2008) Diniz-Filho, J. A. F., Rangel, T. F. L. V. B., & Bini, L. M. (2008). Model selection and information theory in geographical ecology. Global Ecology and Biogeography, 17, 479–488.
  • Dormann & Kaschner (2010) Dormann, C. & Kaschner, K. (2010). Where’s the sperm whale? A species distribution example analysis. http://www.mced-ecology.org/?page_id=355.
  • Dormann et al. (2010) Dormann, C. F., Gruber, B., Winter, M., & Herrmann, D. (2010). Evolution of climate niches in european mammals? Biology Letters, 6, 229–232.
  • Draper (1995) Draper, D. (1995). Assessment and propagation of model uncertainty. Journal of the Royal Statistical Society Series B, 57, 45–97.
  • Elder (2003) Elder, J. F. (2003). The generalization paradox of ensembles. Journal of Computational and Graphical Statistics, 12(4), 853–864.
  • Elith et al. (2006) Elith, J., Graham, C. H., Anderson, R. P., Dudík, M., Ferrier, S., Guisan, A., Hijmans, R. J., Huettmann, F., Leathwick, J. R., Lehmann, A., Li, J., Lohmann, L. G., Loiselle, B. A., Manion, G., Moritz, C., Nakamura, M., Nakazawa, Y., Overton, J. M., Peterson, A. T., Phillips, S. J., Richardson, K., Scachetti-Pereira, R., Schapire, R. E., Soberón, J., Williams, S., Wisz, M. S., Zimmermann, N. E., & Araujo, M. (2006). Novel methods improve prediction of species’ distributions from occurrence data. Ecography, 29(2), 129–151.
  • Elith & Leathwick (2009) Elith, J. & Leathwick, J. R. (2009). Species distribution models: ecological explanation and prediction across space and time. Annual Review of Ecology, Evolution, and Systematics, 40(1), 677–697.
  • Hastie et al. (2009) Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning: Data Mining, Inference, and Prediction, volume 2. Berlin: Springer.
  • Hawkins (2004) Hawkins, D. M. (2004). The problem of overfitting. Journal of Chemical Information and Computer Sciences, 44(1), 1–12.
  • Hegyi & Garamszegi (2011) Hegyi, G. & Garamszegi, L. Z. (2011). Using information theory as a substitute for stepwise regression in ecology and behavior. Behavioral Ecology and Sociobiology, 65(1), 69–76.
  • Horne & Garton (2006) Horne, J. S. & Garton, E. O. (2006). Likelihood cross-validation versus least squares cross-validation for choosing the smoothing parameter in kernel home-range analysis. Journal of Wildlife Management, 70, 641–648.
  • Hurvich & Tsai (1989) Hurvich, C. M. & Tsai, C.-L. (1989). Regression and time series model selection in small samples. Biometrika, 76(2), 297–307.
  • Kullback & Leibler (1951) Kullback, S. & Leibler, R. A. (1951). On information and sufficiency. The Annals of Mathematical Statistics, 1, 79–86.
  • Liaw & Wiener (2002) Liaw, A. & Wiener, M. (2002). Classification and regression by randomForest. R News, 2(3), 18–22.
  • Madigan et al. (1995) Madigan, D., York, J., & Allard, D. (1995). Bayesian graphical models for discrete data. International Statistical Review / Revue Internationale de Statistique, 63(2), 215–232.
  • Mundry (2011) Mundry, R. (2011). Issues in information theory-based statistical inference – commentary from a frequentist’s perspective. Behavioral Ecology and Sociobiology, 65(1), 57–68.
  • Olden et al. (2008) Olden, J. D., Lawler, J. J., & Poff, N. L. (2008). Machine learning methods without tears: A primer for ecologists. The Quarterly Review of Biology, 83(2), 171–193.
  • Raftery (1996) Raftery, A. E. (1996).

    Approximate bayes factors and accounting for model uncertainty in generalised linear models.

    Biometrika, 83(2), 251–266.
  • Recknagel (2001) Recknagel, F. (2001). Applications of machine learning to ecological modelling. Ecological Modelling, 146(1–3), 303–310.
  • Ridgeway et al. (2013) Ridgeway, G. et al. (2013). gbm: Generalized Boosted Regression Models. R package version 2.1.
  • Shen & Huang (2006) Shen, X. & Huang, H.-C. (2006). Optimal model assessment, selection, and combination. Journal of the American Statistical Association, 101(474), 554–568.
  • Shen et al. (2004) Shen, X., Huang, H.-C., & Ye, J. (2004). Adaptive model selection and assessment for exponential family distributions. Technometrics, 46(3), 306–317.
  • Spiegelhalter et al. (2002) Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & van der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society B, 64, 583–639.
  • Stone (1977) Stone, M. (1977). An asymptotic equivalence of choice of model by cross-validation and akaike’s criterion. Journal of the Royal Statistical Society. Series B (Methodological), 39(1), 44–47.
  • Sugiura (1978) Sugiura, N. (1978). Further analysts of the data by akaike’ s information criterion and the finite corrections. Communications in Statistics - Theory and Methods, 7(1), 13–26.
  • Team (2014) Team, R. C. (2014). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
  • Turkheimer et al. (2003) Turkheimer, F. E., Hinz, R., & Cunningham, V. J. (2003). On the undecidability among kinetic models: From model selection to model averaging. Journal of Cerebral Blood Flow & Metabolism, 23(4), 490–498.
  • Venables & Ripley (2002) Venables, W. N. & Ripley, B. D. (2002). Modern Applied Statistics with S. New York: Springer, 4th edition.
  • Wood (2006) Wood, S. (2006). Generalized Additive Models: An Introduction with R. New York: Chapman and Hall/CRC.
  • Wood (2015) Wood, S. N. (2015). Core Statistics. Cambridge: Cambridge University Press.
  • Ye (1998) Ye, J. (1998). On measuring and correcting the effects of data mining and model selection. Journal of the American Statistical Association, 93(441), 120–131.
  • Zhang et al. (2012) Zhang, B., Shen, X., & Mumford, S. L. (2012). Generalized degrees of freedom and adaptive model selection in linear mixed-effects models. Computational Statistics and Data Analysis, 56(3), 574–586.