2. The Setup of the Problem

In model selection, one is concerned with the comparison of families of hypotheses (hereafter “models”). We will reserve the term “hypothesis” for individual members of models. For example, “y = 2x+3+N(0,1)” is a hypothesis belonging to the model {y=ax+b+N(0,1)}. Typically, one’s background theory specifies a finite set of candidate models prior to consulting the data and one is interested in comparing them based on the data. For example, background theory may restrict the set of plausible models to polynomials of degrees no more than 5. One’s goal can then be to find some weighted ordering of such models (here, degrees of polynomials) with respect to various desirable features, such as predictive accuracy or posterior probability.

Suppose you have a suitably large set of observational data consisting of ordered pairs of values of two variables, X and Y, generated by an unknown “true” function, T. The data might, for example, record the length of a metallic bar in different temperatures. Your background theory tells you that Xs and Ys are linearly related (though because of the existence of error, the observed values might not exactly fit a line). You want to find the particular linear function that best fits the data. This is a rather straightforward statistical problem and doesn’t involve model selection, because you are considering only one model (you know X and Y are linearly related). The standard solution on which there is consensus among statisticians is to find the line with maximum likelihood relative to the data, where the likelihood of a hypothesis is defined as the probability (or probability density) of the data conditional on the hypothesis.

A much more difficult question arises if the inference problem concerns the choice between models, say between parabolic functions and linear functions. Here, one cannot simply maximize likelihood, because to do so would almost always lead one to choose a parabolic function (or generally, a member of the largest model). More complicated models have more freedom to fit the data. Thus, the likeliest members of those models are more likely to fit the noise, as opposed to the main pattern, in the data. This is called “over-fitting.” There are various model selection techniques for how to avoid over-fitting. Most of them recommend balancing considerations of goodness-of-fit with data against considerations of simplicity, though the “optimal” balance is naturally different for different techniques, since they either pursue different goals or make different assumptions about the inference problem or both.

Before we proceed, some notational conventions must be mentioned. In order to avoid confusion, we refer to random variables by capital letters, single pieces of data by lower-case letters and data sets by bolded lower-case letters. In order to distinguish between specific data sets and data sets considered as random variables, we show the former by subscripted, bolded lower-case letters, like y ₀ , and the latter by non-subscripted, bolded lower-case letters, like, y. Consider a model, F(θ), defined over a parameter space Θ. For example, if F is the family of linear functions with a normally-distributed error with mean 0 and variance 1, it can be characterized as F: {y = ax+b+N(0,1); (a,b) ∈ R²}. Here the parameter space of F is R².Footnote ³ We can associate a likelihood function, ℒ(θ) = ℒ(θ,z) = g(z|θ), to F, where g(z |θ) is the probability density of obtaining the n-tuple data set z ={z_1,z₂,…,z_n} [where z_i is the ith observed datum (x _i ,y _i )] conditional on θ being the true parametric value. The value of θ for which ℒ(θ) is maximized is called the Maximum Likelihood Estimate (MLE) of F, and we denote it by $\hat \theta $ . We denote the member of F obtained by taking θ = $\hat \theta $ by L(F).Footnote ⁴ Note that both $\hat \theta $ and L(F) are functions of data; that is, $\hat \theta $ = $\hat \theta $ (y) and L(F) = L(F,y), where y is the data.

In the Akaikean framework, the goal is to find a plausible estimate of the predictive accuracies of various models. The predictive accuracy of a model, M, is a measure of how close, on average, the predictions of the best fitting member of M, with respect to an initial data set, are to subsequent data generated by the same generating function. Suppose you obtain a data set, y, generated by T and use it to determine L(M,y). Now obtain a new data set, x, and calculate the logarithm of its likelihood (hereafter, “log-likelihood”) of L(M,y) with respect to x. The average value—with respect to both x and y—of this log-likelihood is called the predictive accuracy of M. The following defines A(M), the predictive accuracy of model M.Footnote ⁵

(1) $${\rm{A}}\left( M \right){ = ^{\rm df}}{\rm{E_y}{E_x}[\log (g(x|\hat \theta \left( y \right)))]}{\rm{ }} = {E_{\hat \theta }}{\rm{E_x}[\log (g(x|\hat \theta ))]}$$

E _y (.) and E _x (.) are both expectations with respect to T. The last term on the right-hand side better captures the nature of predictive accuracy, by emphasizing that it averages over the MLE.Footnote ⁶

Akaike showed that the AIC value of a model can be used to provide an estimate of its predictive accuracy. AIC of model F is defined thus.

(2) $${\rm{AIC}}\left( F \right){ = ^{{\rm{df}}}} - 2{\rm{log{\scr L}}}({\rm{L}}\left( F \right)) + 2k$$

logℒ( ${\rm{L}}\left( F \right)$ ) is the log-likelihood of L(F). If error is normally distributed, this becomes the familiar sum of squares of error terms of L(F). k is the dimension of the parameter space of F, which is usually equal to its number of adjustable parameters.

If one is interested in finding the model with the highest predictive accuracy, the Akaikean frameworkFootnote ⁷ recommends choosing the model with minimum AIC. However, as Forster and Sober observe, blindly minimizing AIC is problematic. Suppose one wishes to find the predictively most accurate hypothesis among polynomials of degrees 3 or less. One must find the family with minimum AIC and then select the likeliest member of that family. Let M _i be the family of polynomials of degree i-1 with i adjustable parameters. Suppose M ₃ turns out to have the lowest AIC value. Since the family of parabolic functions is embedded in the family of cubic functions, the likeliest member of the cubic family (L(M ₄ )) has a better (or equally good) fit with the data than the likeliest member of the parabolic family (L(M ₃ )). Now construct an ad hoc family, {L(M ₄ )}, whose only member is L(M ₄ ). The number of adjustable parameters in {L(M ₄ )} is 0 (because it is a singleton family) and therefore, its AIC value is equal to -2logℒ(L(M ₄ )). We have, AIC({L(M ₄ )}) = -2logℒ(L(M ₄ )) $ \le $ -2logℒ(L(M ₃ )) < AIC(M ₃ ), where logℒ(L(M ₄ )) is the log-likelihood of the best fitting member of M₄. Indeed, {L(M ₄ )} has the lowest possible AIC value (relative to the extant data) among all families the members of which are restricted to polynomials of degrees 3 or less. Thus, if we blindly minimize AIC scores, we must choose {L(M ₄ )} (hereafter, “the sub-family model”) as the predictively most accurate model, which is tantamount to choosing the likeliest hypothesis at our disposal and giving no weight to simplicity. This is what Forster and Sober have called the sub-family problem.Footnote ⁸

In order to show the importance of IMR, we will contrast the sub-family problem with another inference problem, which is very similar but differs in only one salient way: in that problem, IMR is respected. Suppose you have reasons independently of the extant data (e.g., theoretic reasons) for including the singleton model {5x³+6x²+2x+7+N(0,1)} among your candidate models. Then you obtain a data set, y ₀ , and you observe that L(M ₄ ,y ₀ )=5x³+6x²+2x+7+N(0,1).Footnote ⁹ This inference problem, which we will call the singleton family problem,Footnote ¹⁰ is very similar to the sub-family problem. The set of candidate models (M ₁ ,M ₂ ,M ₃ ,M _4, {L(M ₄ ,y ₀ )}) and the data (y ₀ ) are the same. However, they have an important difference: here IMR is respected because you had reasons to include {5x³+6x²+2x+7+N(0,1)} (which happens to be identical with {L(M ₄ ,y ₀ )}) among your candidate models. As we shall see, this makes a big difference.

3. BIC

In the BIC-based framework, the average likelihood of a modelFootnote ¹¹ is estimated by −½ times exponential of its BIC defined as follows,

(3) $${\rm{BIC}}\left( F \right){ = ^{{\rm{df}}}} - 2{\rm{log{\scr L}}}({\rm{L}}\left( F \right)) + k{\rm{log}}\left( n \right),$$

where n is the number of data points. Lower BIC scores are better. Since this is a Bayesian framework, its ultimate goal is to determine the posterior probabilities of candidate models. It is customary—though by no means necessary—to assign equal prior probabilities to all candidate models. Thus, the most probable model is often the one with the lowest BIC value. Now, the sub-family model has the lowest possible BIC value, since it has no adjustable parameters and its only member has maximum log-likelihood. Thus, we are faced with the sub-family problem: the sub-family model appears to be the most probable model.

A referee suggests that some of the general problems for the BIC-based framework might complicate our solution to the sub-family problem. Thus, before offering our solution, we will briefly mention one of those problems, which is most relevant. If the models are nested (as in our example of polynomial models), then larger models entail smaller ones. (The set of polynomials of degree n contains the set of polynomials of degree n-1.) It follows that P(M _n ) ≥ P(M _n-1 ), no matter what the data is. If so, it is difficult, in this framework, to make sense of the fact that scientists sometimes prefer smaller models to larger ones. To the best of our knowledge, there has been no fully satisfactory response to this problem. Forster and Sober (Reference Forster and Sober1994) discuss the following way to address this difficulty. Instead of M ₂ , construct M ₂ * = M ₂ - M ₁ (and so on for larger models), so that no model entails any other. Then compare those newly-formulated models. Forster and Sober find this maneuver unsatisfactory, because it changes the subject. The question was why scientists prefer M ₁ to M ₂ , not M ₁ to M ₂ *. This is a fair objection, a proper reply to which (if it exists) goes beyond the scope of this essay. Here we use this mathematical maneuver to offer a solution for the sub-family problem, but we don’t claim to offer a solution for the above problem or any other general problem for BIC. All we wish to establish is that if BIC can be rescued from the general difficulties it faces, it won’t be further subject to the sub-family problem.

Construct the non-nested models (M ₁ -M ₄ *) as described above. Make sure that the prior probability functions over the members of your non-nested models do not have probability masses. That is, if π₄( ${\rm{\theta }}$ ) is the prior probability function over the members of M ₄ *, then make sure for no single value of ${\rm{\theta }}$ , π₄( ${\rm{\theta }}$ ) > 0. (If this is not the case, BIC is not a good approximation of average likelihood. This is another limitation of BIC, again, independently of the sub-family problem.) How can you make sure this is the case? If for ${{\rm{\theta }}_1}$ , π₄( ${{\rm{\theta }}_1}$ ) > 0, construct the singleton model { ${{\rm{\theta }}_1}$ }, where P({ ${{\rm{\theta }}_1}$ }) = π₄( ${{\rm{\theta }}_1}$ ). Then redefine M ₄ * in the following way, M ₄ * _new = M ₄ * _old -{ ${{\rm{\theta }}_1}$ }, with P(M ₄ * _new ) = P(M ₄ * _old ) - P({ ${{\rm{\theta }}_1}$ }). Once this is done, BIC can be used as an estimate of average likelihood of the non-nested models. Now suppose we obtain a data set y ₀ and {L(M ₄ ,y ₀ )} = {5x³+6x²+2x+7+N(0,1)}. There are two possibilities. Either P({5x³+6x²+2x+7+N(0,1)}) > 0, or P({5x³+6x²+2x+7+N(0,1)}) = 0, where P(.) is the prior probability function.

If P({5x³+6x²+2x+7+N(0,1)}) > 0, then you had a positive reason to include {5x³+6x²+2x+7+N(0,1)} among your candidate models independently of the dataFootnote ¹² and this is an instance of the singleton family problem—considering ({5x³+6x²+2x+7+N(0,1)}) does not violate IMR. The likelihood of a model is proportional to the exponential of −½BIC and since {5x³+6x²+2x+7+N(0,1)} has an exceptionally low BIC, its likelihood will be massively higher than other models. Thus, unless its prior probability is extremely lower than other models, it will be by far the most probable model among the non-nested ones. This is a welcome result. If one had theoretic reasons that a single parameter value has positive probability, and one subsequently learns that this single hypothesis fits data excellently well, the data provides very powerful evidence for that hypothesis and ought to make one significantly more confident of its truth. Again, this doesn’t solve the above-mentioned general difficulty about BIC. In the singleton family problem, {5x³+6x²+2x+7+N(0,1)} might end up having a higher posterior probability than M ₄ ^* but it can never have a higher probability than M ₄ .

The sub-family problem corresponds to the case where P({5x³+6x²+2x+7+N(0,1)}) = 0. Here you have no reason to include {5x³+6x²+2x+7+N(0,1)} among the candidate models independently of the data. Thus, to choose {5x³+6x²+2x+7+N(0,1)} because of its excellent BIC score involves violating IMR. The BIC-based framework disallows choosing {5x³+6x²+2x+7+N(0,1)}, because if P({5x³+6x²+2x+7+N(0,1)})=0, no matter how likely {5x³+6x²+2x+7+N(0,1)} is (no matter how good its BIC score is), its posterior probability will remain zero. As we shall see, although the Akaikean framework doesn’t appeal to model priors, it equally presupposes IMR.

4. Forster and Sober’s Solution

For large data sets, AIC is an approximately unbiased estimator of predictive accuracy. An estimator is statistically unbiased if its expected value equals the value it estimates. The following equation expresses this fact.

(4) $${{\rm{E}}_{\bf{y}}}\left[ { - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} {\rm{AIC}}\left( {M,{\bf{y}}} \right)} \right] = {\rm{A}}\left( M \right),$$

where y is a random data set. In this paper, we are not concerned with the approximate nature of this equation. So we will talk about the “unbiasedness” of AIC for convenience.

Forster and Sober argue that the AIC of the sub-family model is an unbiased estimator of its predictive accuracy, but statistical unbiasedness is not the only criterion by which to judge an estimate. They introduce another such criterion called “epistemic unbiasedness” by the following example. Consider a simple measurement of an object’s mass with a kitchen scale. Normally, the measured value is an unbiased estimator of the actual value, because it is just as likely to over-measure the mass by a given amount as it is to under-measure it by that same amount.

It is helpful to make a distinction between an “estimator” and an “estimate” here. An estimator is a function of data that yields individual estimates. Statistical unbiasedness is a feature of an estimator, while epistemic bias is a feature of an individual estimate. Forster and Sober argue that AIC of the sub-family model is statistically unbiased (qua estimator) but epistemically biased (qua estimate).Footnote ¹³ They argue that the AIC of the sub-family model is not a good estimate because it is epistemically biased. To show this, they appeal to the following “meta-theorem” about AIC.

Error[Estimated(A(F))] =^df A(F) – (−½)AIC(F) = Residual Fitting Error + Common Error + Sub-family Error (Forster and Sober Reference Forster and Sober1994, 19).

This theorem concerns the error involved in taking −½AIC as an estimate of predictive accuracy. Forster and Sober argue that the first two terms on the right-hand side are both statistically and epistemically unbiased, but the third term, although statistically unbiased, is sometimes epistemically biased, which (given the epistemic unbiasedness of the other two terms) makes the total error sometimes epistemically biased. They further argue that an important occasion in which this happens is in the sub-family problem. Thus, a fuller understanding of the Akaikean framework, which includes this meta-theorem, dissolves the sub-family problem.

Here is why Forster and Sober believe the sub-family error is epistemically biased for the sub-family model. Suppose we embed the parameter spaces of all our models in a larger parameter space (call it K) that contains the truth, T. Forster and Sober state that this space can be considered as a vector space in such a way that i) the closer a point in this space is to truth, the higher its predictive accuracy; and ii) the sub-family error is equal to the scalar product of the following two vectors in this space: the vector that goes from T to the likeliest hypothesis in K, L(K), (i.e., $\overrightarrow {T.{\rm{L}}\left( K \right)} $ ) and the vector $\overrightarrow {T.{{\rm{\theta }}_0}} $ that goes from T to the (unknown) predictively most accurate member of M, which we denote by θ₀.Footnote ¹⁴ These vectors are shown in figure 1 below.Footnote ¹⁵

Figure 1.

T is the truth. L(K) is the likeliest hypothesis in K, which contains T. For a Gaussian distribution of error, L(K) is equally likely to fall on any point on the circle. T and L(K) are independent of the model. θ₀ is the predictively most accurate hypothesis in the model (i.e., the closest member of the model to T).

All three points are in general unknown. However, the scalar product of $\overrightarrow {T.{\rm{L}}\left( K \right)} $ and $\overrightarrow {T.{{\rm{\theta }}_0}} $ is equal to the product of their lengths multiplied by the cosine of the angle between them. Forster and Sober argue that for the sub-family model, the tips of the two vectors tend to be close. Here is their argument.

On Forster and Sober’s view, we shouldn’t follow the sub-family policy, because if we do, the AIC values of the models we construct, although statistically unbiased, tend to be epistemically biased (over)estimates of the predictive accuracies of these models.

Before discussing what we find potentially problematic about this argument, we would like to give a summary of what we will claim, in order to avoid confusion. When Forster and Sober say that the AIC of the sub-family model is statistically unbiased, they are referring to the fact that the estimate AIC({L(M ₄ ,y ₀ )},y ₀ ) is an instantiation of the estimator AIC({L(M ₄ ,y ₀ )},y), which is—we agree—an unbiased estimator. We also agree that in the sub-family problem, AIC({L(M ₄ ,y ₀ )},y ₀ ) is epistemically biased. However, we suggest that a better estimator (than AIC({L(M ₄ ,y ₀ )},y)) for determining the (de)merits of the estimate AIC({L(M ₄ ,y ₀ )},y ₀ ) in the sub-family problem is AIC({L(M ₄ ,y)},y), which is statistically biased. This is because in the sub-family problem, AIC({L(M ₄ ,y ₀ )},y ₀ ) is more similar to other instantiations of AIC({L(M ₄ ,y)},y) than to other instantiations of AIC({L(M ₄ ,y ₀ )},y). The exact opposite situation holds for the singleton family problem. There the statistical bias of AIC({L(M ₄ ,y ₀ )},y) provides more relevant information (than that of AIC({L(M ₄ ,y)},y)) for deciding how good the estimate AIC({L(M ₄ ,y ₀ )},y ₀ ) is. Therefore, we will argue that the Akaikean framework treats the sub-family problem and the singleton family problem differently.

Now we will unpack this. Consider the singleton family problem first. That is, suppose we had reasons independently of y ₀ for considering the model {5x³+6x²+2x+7+N(0,1)}. Then we obtained y ₀ and observed that L(M ₄ ,y ₀ )=5x³+6x²+2x+7+N(0,1). Thus, including {5x³+6x²+2x+7+N(0,1)} among the models doesn’t violate IMR. Is AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) an epistemically biased estimate in this case, according to Forster and Sober? We don’t know for sure because they don’t discuss this problem. However, if we are to infer from the literal meaning of their argument, they would have to say “yes,” because in every relevant respect to their solution, the two problems are identical. In their argument, the relevant factors are the low-dimensionality of the model and its closeness to data, which are shared in the two problems. If they are committed to this idea, then this is the only point of disagreement between our account and theirs. We believe that in the singleton family problem, the exceptionally low AIC score of {5x³+6x²+2x+7+N(0,1)} is exceptionally good reason that {5x³+6x²+2x+7+N(0,1)} has high predictive accuracy. There is nothing “unfair” in the AIC score of a low-dimensional model with good fit per se. The AIC scores of such models are “too good to be true” only if the model has these features because it was designed in an ad hoc fashion to have a low AIC score.

For this to be true, something must be missing in Forster and Sober’s account of the sub-family error as applied to the singleton family problem. Here is what is missing. They offer a consideration (hereafter consideration₁) that ${{\rm{\theta }}_0}$ tends to be close to L(K), which leads to the sub-family error for {5x³+6x²+2x+7+N(0,1)} being large and positive. However, there is a competing consideration (hereafter, consideration₂) that mitigates the effect of consideration₁: for models with low AIC scores, ${{\rm{\theta }}_0}$ tends to be close to T. (A low AIC score means high predictive accuracy, which means closeness to truth.) Consideration₂ is a reason for sub-family error to be small, because the error is equal to the product of the lengths of the two vectors $\overrightarrow {T.{\rm{L}}\left( K \right)} $ and $\overrightarrow {T.{{\rm{\theta }}_0}} $ times the cosine of the angle between them. Consideration₂ is a reason that the length of $\overrightarrow {T.{{\rm{\theta }}_0}} $ is small. In both the sub-family problem and the singleton family problem, AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) is exceptionally low, which is a reason to think that 5x³+6x²+2x+7+N(0,1) is close to truth. This is the case, unless one otherwise knows that AIC({5x³+6x²+2x+7+N(0,1)}) is not a good estimate of the predictive accuracy of {5x³+6x²+2x+7+N(0,1)}. We will argue in the next section that in the sub-family problem, one knows this independently. Thus, consideration₂ is irrelevant to the sub-family problem and AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) is epistemically biased in that problem. However, our argument in the next section doesn’t apply to the singleton family problem. In that problem, we are left with two competing considerations bearing on how good an estimate AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) is, and in general we have no way of comparing the relative strengths of these considerations. Notice that consideration₂ doesn’t tell us anything about the sign of the sub-family error (whether it is negative or positive). Thus, we must expect to have a positive epistemic error but smaller in absolute value relative to the sub-family problem.

The amount of bias matters a lot. In both the sub-family problem and the singleton family problem, AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) is significantly better than the AIC of other models. For example, AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) =AIC(M ₄ ,y ₀ ) – 4,Footnote ¹⁷ and 4 units of AIC difference is usually a big difference. (To have an intuitive idea for why this is the case, consider the fact that if the data is perfectly linear, M ₂ will have only a 2 units AIC advantage relative to M ₄ .) In the sub-family problem, the average sub-family error is exactly 4, which means that the entire difference in AIC score (between ({5x³+6x²+2x+7+N(0,1)} and M ₄ ) is due to error. However, in the singleton family problem, little of the AIC difference between {5x³+6x²+2x+7+N(0,1)} and M ₄ (which is 4 units) is due to epistemic error. Therefore, {5x³+6x²+2x+7+N(0,1)} is in great shape relative to M ₄ . This is a welcome result. If you have a singleton model among the candidate models for reasons independently of the extant data, and it fits the data as precisely as the family of cubic functions, you ought to be very confident of its predictive accuracy even if its AIC score is slightly biased.

The following numerical example can help illustrate our point. Here all the data is generated by the function T= 0.2x⁵–0.2x⁴–3x³+x²–1+N(0,1). Figure 2a depicts the sub-family error of the model {L(M ₄ ,y)} for 10⁷ values of y each consisting of 100 data points.Footnote ¹⁸ That is, for each data set, y, L(M ₄ ,y) is determined separately and then its corresponding sub-family error is calculated.Footnote ¹⁹ This is essentially a 10⁷ repetition of the sub-family problem. Evidently, the sub-family error of {L(M ₄ ,y)} tends to be positive and large, as Forster and Sober rightly argue. Figure 2b depicts the sub-family error of the fixed singleton model M ₀ :{–3x³+x²–1+N(0,1)}, which is chosen because it is close to T. The sub-family error of this fixed model is statistically unbiased, again as observed by Forster and Sober.

Figure 2.

Histogram of the Sub-family Error of {L(M ₄ ,y)}. Histogram of the Sub-family Error for M ₀ :{–3x³+x²–1+N(0,1)}.

The difficult part is how to simulate the singleton family problem, because in that problem you must have a singleton family for reasons independently of the data and then something amazing happens: the only member of that model turns out to be equal to the best fitting member of your largest model. (Even if your singleton family contains the truth, this is unlikely to happen, unless the data is huge.) However, we can approximate this situation. Figure 3a is the histogram of the distance between M ₀ and {L(M ₄ ,y)} for the 10⁷ data sets. Instead of looking at cases where {–3x³+x²–1+N(0,1)}={L(M ₄ ,y)}, we first looked at cases in which M ₀ is “close” to {L(M ₄ ,y)}. This corresponds to cases where the member of the singleton family is not exactly identical with L(M ₄ ,y) but is close to it (both in the parameter space and in terms of log-likelihood). In 3b, the sub-family error for M ₀ is depicted only for those cases in which ||M ₀ -L(M ₄ ,y)||< 0.3 (cases that are left of the red line in 3a). The choice of 0.3 is arbitrary. We chose this value so that we can still have a large number of cases. (Smaller values result in larger errors because for them consideration₁ becomes stronger, but due to the small number of cases, the histograms become very jagged.)

Figure 3.

Histogram of the distance between M ₀ and L(M ₄ ,y). Histogram of the sub-family error for M ₀ only for those data sets for which ||M ₀ -L(M ₄ ,y)||< 0.3 (left of the red line in 3a).

The average sub-family error for cases depicted in figure 3b is 0.4286. When we took this condition to the limit (i.e., ||M ₀ -L(M ₄ ,y)||→0), the error approached a value slightly less than 1. A comparison between 2a and 3b shows the mitigating effect of consideration₂. In 2a, the average L(M ₄ ,y) is not particularly close to T; thus, the only relevant consideration is consideration₁. The average sub-family error was 3.9991 in 2a.Footnote ²⁰ In 3b, the model is still singleton and close to the data (thus consideration₁ is still pertinent), but because M ₀ is close to T, the average sub-family error was 0.4286, quite smaller than 3.9991. The importance of this fact can be best understood with the help of figure 4, which depicts the histogram of AIC(M ₄ )-AIC(M ₀ ) for those cases depicted in 3b. For ||M ₀ -L(M ₄ ,y)||< 0.3, the average AIC difference was 3.2197.

Figure 4.

Histogram of the AIC difference between M ₀ and M ₄ for data sets for which ||M ₀ -L(M ₄ ,y)||< 0.3 (left of the red line in 3a).

The important point is that as we focus on data sets for which M ₀ is closer and closer to L(M ₄ ,y), AIC(M ₀ ) becomes smaller and smaller (better) but at the same time the sub-family error becomes larger and larger. For ||M ₀ -L(M ₄ ,y)||→0, (AIC(M ₀ ,y)–AIC(M ₄ ,y))→ -4 and the sub-family error approached 1 in this particular case. That is, if you correct for the sub-family error, AIC(M ₀ ,y) is still 3 units better than that of AIC(M ₄ ,y). For the average case that satisfied ||M ₀ -L(M ₄ ,y)||< 0.3, if you correct for the sub-family error (i.e., subtract the average sub-family error, 0.4286), you still have a better AIC value of 2.7911. This illustrates the main point of our argument. In the singleton family problem, our model is singleton and very close to the data. However, its excellent AIC score (although slightly biased) is excellent evidence that it is predictively accurate.

Objection: the difference between the sub-family problem and the singleton family problem is only a historical fact about how the model was constructed. So it cannot affect how good an estimate AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) is.

Answer: certain facts about AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) are not affected by this historical fact, including the definition of A({5x³+6x²+2x+7+N(0,1)}), the value of AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) itself and the statistical bias of the estimator AIC({5x³+6x²+2x+7+N(0,1)},y). However, this list does not exhaust all the relevant information bearing on how good an estimate AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) is. Indeed, Forster and Sober introduced the notion of epistemic bias in order to be able to account for the intuitive idea that in the sub-family problem there is something wrong with AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) as an estimate, which cannot be captured by the items in the above list. We think what is wrong with this estimate is that in the sub-family problem, our model is tailored to y ₀ Footnote ²¹ in an ad hoc fashion in order to have an optimal AIC score. However, this is not the case in the singleton family problem. In that problem, {5x³+6x²+2x+7+N(0,1)} is a contender model for independent reasons. And that such a simple model fits y ₀ so well (without us having chosen it because it fits y ₀ well) is excellent reason that it is highly predictively accurate. In other words, the historical fact that constitutes the difference between the two problems bears important information about the ad hocness of AIC({5x ³ +6x ² +2x+7+N(0,1)},y ₀ ) as an estimate. Of course, this intuitive idea needs a technical explanation. That is what we will offer in the next section.

5. Our Solution

We think a better solution for the sub-family problem can be given by studying another estimator namely AIC({L(M ₄ ,y)},y) (instead of AIC({L(M ₄ ,y ₀ )},y)) one of the instantiations of which is the estimate AIC({L(M ₄ ,y ₀ )},y ₀ ). But how can one decide which estimator provides better information about the merits/demerits of the individual estimate? In order to answer this question, consider why statisticians study the statistical biases of estimators in the first place. Suppose a ₁ is an estimate of the quantity, a. For example, a ₁ might be the value one reads on a kitchen scale when one weighs an apple. Naturally, one might wish to know how good an estimate a ₁ is, but usually one cannot directly talk about how good or bad a single estimate is. However, sometimes one can talk about how good, in general, the estimates obtained by the “same” procedure tend to be. That is, one can talk about various features of an estimator by considering repeated estimates obtained by the “same” procedure. Statistical bias is one such feature. An estimator whose expected value equals the value it estimates is statistically unbiased. But what exactly does this tell us about the individual estimate? Insofar as the individual estimate is similar to the other instantiations of the estimator, the expected error of the estimator (its statistical bias) contains information about how good the estimate is. Importantly, the inherent vagueness in what it means for the procedure to be the “same” makes it the case that a particular estimate can be an instantiation of more than one estimator.Footnote ²² However, it doesn’t follow that the statistical biases of those estimators bear equally valuable information on the merits of the individual estimate. If the estimate a ₁ is an instantiation of two estimators A and A*, and if other instantiations of A better resemble a ₁ than other instantiations of A* in ways that affect the value of the estimate, then the statistical bias of A bears more pertinent information than the statistical bias of A* on the merits/demerits of a ₁ .

In ordinary inference problems, there is only one estimator one would naturally associate with an individual AIC score (qua estimate). Thus, these considerations are usually unimportant. However in the sub-family problem, things are different because not only the individual AIC score is a function of the data, but the model itself is designed on the basis of the data too. Thus, the AIC score one calculates - AIC({L(M ₄ ,y ₀ )},y ₀ ) - is an instantiation of two estimators: AIC({L(M ₄ ,y ₀ )},y) and AIC({L(M ₄ ,y)},y). Notice that both of these estimators can be understood as the AIC of “the sub-family model.” If you understand “the sub-family model” rigidly to refer to the model constructed after obtaining y ₀ , you’ll have AIC({L(M ₄ ,y ₀ )},y); and if you understand it nonrigidly to refer to whichever model one constructs on the basis of data, you’ll have AIC({L(M ₄ ,y)},y). The former is an unbiased estimator, as Forster and Sober correctly claim. The latter is a biased estimator, as we shall argue later. But which one gives us better information about the merits/demerits of the individual estimate AIC({L(M ₄ ,y ₀ )},y ₀ )? In order to answer this question, consider another instantiation of each estimator, say for data set y ₁ . The estimator AIC({L(M ₄ ,y ₀ )},y) yields AIC({L(M ₄ ,y ₀ )},y ₁ ). This keeps the model {5x³+6x²+2x+7+N(0,1)} and calculates the AIC score of that fixed model with respect to y ₁ . However AIC({L(M ₄ ,y)},y) yields AIC({L(M ₄ ,y ₁ )},y ₁ ). Here one constructs the model {L(M ₄ ,y ₁ )} and calculates AIC({L(M ₄ ,y ₁ )}) with respect to y ₁ . AIC({L(M ₄ ,y ₀ )},y ₁ ) is dissimilar to AIC({L(M ₄ ,y ₀ )},y ₀ ) in the sub-family problem in an important respect. In AIC({L(M ₄ ,y ₀ )},y ₁ ) the model is chosen because it fits one data set (y ₀ ) well but its AIC score is calculated with respect to another data set (y ₁ ). Whereas in AIC({L(M ₄ ,y ₀ )},y ₀ ), the model is chosen because it has excellent fit with y ₀ and its AIC score is calculated with respect to that same data set. Obviously, AIC({L(M ₄ ,y ₀ )},y ₀ ) and AIC({L(M ₄ ,y ₁ )},y ₁ ) are similar in this respect. Another way one can see this point is by considering the following question: What would have been the AIC estimate if instead of y ₀ one had obtained y ₁ in the sub-family problem? Clearly the answer is AIC({L(M ₄ ,y ₁ )},y ₁ ). This is the reason we believe the statistical bias of the estimator AIC({L(M ₄ ,y)},y) provides better information about the estimate AIC({L(M ₄ ,y ₀ )},y ₀ ) in the sub-family problem.Footnote ²³

What about the singleton family problem? Things are quite different in that problem. What would have been the AIC estimate if instead of y ₀ one had obtained y ₁ in the singleton family problem? Clearly AIC({5x³+6x²+2x+7+N(0,1)},y ₁ ), which happens to be equal to AIC({L(M ₄ ,y ₀ )},y ₁ ), because here {5x³+6x²+2x+7+N(0,1)} is a contender model regardless of the fact that 5x³+6x²+2x+7+N(0,1) = L(M ₄ ,y ₀ ). Therefore, for this problem the estimator AIC({L(M ₄ ,y)},y) is irrelevant, because we care about the singleton family {5x³+6x²+2x+7+N(0,1)}. Here the statistical bias of AIC({L(M ₄ ,y ₀ )},y) gives one the relevant information on how good an estimate AIC({5x³+6x²+2x+7+N(0,1)},y ₀ ) (or incidentally, AIC({L(M ₄ ,y ₀ )},y ₀ )) is.

Now we would like to show that AIC({L(M ₄ ,y)},y) is in fact a biased estimator. This might appear to contradict Akaike’s results, but it doesn’t. Those results presuppose IMR according to which the set of candidate models must be determined independently of the current data. In fact, Akaike’s results about models that respect IMR help us prove that AIC({L(M ₄ ,y)},y) is biased.

In order to show this, we first show that the expectation of −½(AIC({L(M ₄ ,y)},y) with respect to y is larger than the predictive accuracy of any singleton family whose only member is a member of M ₄ . Suppose f is a variable that ranges over the members of M. For each value of f, we can define the singleton family {f}. Also suppose m is an arbitrary member of M and {m} is the singleton family whose only member is m. By definition of L(M ₄ ,y) we have,

(5) $$ - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} {\rm{AIC}}\left( {\left\{ {{\rm{L}}\left( {{M_4},{\bf{y}}} \right)} \right\},{\bf{y}}} \right) = {\rm{max}}[ - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} ({\rm{AIC}}\left( {\left\{ {\rm{f}} \right\},{\bf{y}}} \right),{\rm{f}} \in {M_4}] \ge - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} {\rm{AIC}}\left( {\left\{ {\rm{m}} \right\},{\bf{y}}} \right)$$

If we take expectation with respect to y we have:

(6) $$\matrix{ {{{\rm{E}}_{\bf{y}}}[ - {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}{\rm{AIC}}(\left( {\left\{ {{\rm{L}}\left( {{M_4},{\bf{y}}} \right)} \right\},{\bf{y}}} \right)] = {{\rm{E}}_{\bf{y}}}\left[ {{\rm{max}}} \right[ - {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}{\rm{AIC}}\left( {\left\{ {\rm{f}} \right\},{\bf{y}}} \right),{\rm{f}}} \hfill \cr {\quad \quad \in {M_4}]] \ge {{\rm{E}}_{\bf{y}}}\left[ { - {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}{\rm{AIC}}\left( {\left\{ {\rm{m}} \right\},{\bf{y}}} \right)} \right] = {\rm{A}}\left( {\{ {\rm{m}}\} } \right)} \hfill \cr } $$

The last equality obtains because of (4). In (6), only if m = L(M ₄ ) the equality holds. However, L(M ₄ ) is a function of data. For any m, such that m = L(M ₄ ,y _i ), there is (almost always) a data set y _j generated by the same generating function such that m is not the likeliest member of M ₄ with respect to y _j . Thus, in fact we have a stronger result than (6):

(7) $$\matrix{ {{{\rm{E}}_{\bf{y}}}[ - {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}{\rm{AIC}}\left( {\left\{ {{\rm{L}}\left( {{M_4},{\bf{y}}} \right)} \right\},{\bf{y}}} \right)] = {{\rm{E}}_{\bf{y}}}\left[ {{\rm{max}}} \right[ - {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}{\rm{AIC}}\left( {\left\{ {\rm{f}} \right\},{\bf{y}}} \right),{\rm{f}}} \hfill \cr {\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \in {M_4}\left] {\left] {{\rm{ }}{{\rm{E}}_{\bf{y}}}} \right[ - {\raise0.5ex\hbox{$\scriptstyle 1$} \kern-0.1em/\kern-0.15em \lower0.25ex\hbox{$\scriptstyle 2$}}{\rm{AIC}}\left( {\left\{ {\rm{m}} \right\},{\bf{y}}} \right)} \right] = A\left( {\{ {\rm{m}}\} } \right)} \hfill \cr } $$

It follows from (7) that the average value of −½AIC({L(M ₄ ,y),y}) is strictly larger than the predictive accuracy of any singleton family one can construct from the members of M ₄ (since m is an arbitrary member of M ₄ )—including, of course, {L(M ₄ ,y ₀ )}.

Since {L(M ₄ ,y)} is a random variable and a function of data, A({L(M ₄ ,y)}) is a random variable too and will vary for different data sets. This is unlike the usual application of the Akaikean framework, where the model is fixed and its predictive accuracy is a fixed number. Indeed if {L(M ₄ )} was not a function of y, then by (4), −½AIC({L(M ₄ )}) would have been an unbiased estimator of A({L(M ₄ )}). However, regardless of the data at hand, since L(M ₄ ) is a member of M₄, we have A({L(M ₄ ,y)}) ≤ max[(A({f}), f ∈ M ₄ ]. And since (7) is true for all m, m ∈ M ₄ , then we have,

(8) $${{\rm{E}}_{\bf{y}}}\left[ { - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} {\rm{AIC}}\left( {\left\{ {{\rm{L}}\left( {{M_4},{\bf{y}}} \right)} \right\},{\bf{y}}} \right)} \right] \gt {\rm{A}}\left( {\{ {\rm{L}}\left( {{M_4},{\bf{y}}} \right)\} } \right) $$

That is, AIC({L(M ₄ ,y)},y) is statistically biased. A comparison between equations (4) and (8) shows the difference between AIC({L(M ₄ )}) and AIC of ‘normal’ models that are constructed independently of the data. In the same way that equation (4) motivates using AIC scores of ‘normal’ models as an estimate of their predictive accuracies, equation (8) shows why −½ times AIC({L(M ₄ ,y ₀ )},y ₀ ) is not a good estimate of A({L(M ₄ ,y ₀ )}).

There is another way of looking at AIC({L(M ₄ ,y)},y) that makes it obvious why it is a biased estimator. Since {L(M ₄ ,y)}is singleton,

(9) $${\rm{AIC}}\left( {\left\{ {{\rm{L}}\left( {{M_4},{\bf{y}}} \right)} \right\},{\bf{y}}} \right) = - 2{\rm{log\scr L}}\left( {\left\{ {{\rm{L}}\left( {{M_4},{\bf{y}}} \right)} \right\},{\bf{y}}} \right) = {\rm{AIC}}\left( {{M_4},{\bf{y}}} \right) - 4,$$

The last equality obtains because M ₄ has 4 adjustable parameters. But by the definition of predictive accuracy (equation (1)) for model M ₄ , the average predictive accuracy of {L(M ₄ ,y)} (i.e., E _y [A({L(M ₄ )})]) must equal A(M ₄ ). It follows that AIC({L(M ₄ ,y)},y) is on average 4 units lower than A({L(M ₄ ,y)}). The fact that the estimator AIC({L(M ₄ ,y)},y) is biased gives us information about why the individual estimate AIC({L(M ₄ ,y ₀ )},y ₀ ) is not a good estimate of predictive accuracy. This concludes our technical treatment of the sub-family problem. In the next section, we will talk about another inference problem in which IMR is violated.

7. Too Many Models

Post hoc model construction is not the only way one can violate IMR. Another way is by adding models to one’s set of candidate models without any reason for them to be considered. This simply increases the probability that a model with low predictive accuracy will have a good AIC score because it fits the current data well. An extreme version of this strategy is to add every member of the largest candidate model as a singleton model. In our working example of polynomial models, this would involve adding a singleton family {M_a,b,c,d: y = ax³+bx²+cx+d+N(0,1)} for all (a,b,c,d) ∈ R⁴. The number of candidate models will be infinite, but the one with the lowest AIC score is obviously {L(M ₄ ,y)}. Therefore, in this case one would always choose the same model as one does in the sub-family problem. No doubt this is a problematic practice, but the question is why?

In order to see why, first consider the rather common practice of computing AIC scores for all candidate models, choosing the model with minimum AIC as the winner and taking its AIC value as an estimate of its predictive accuracy. Call this policy the minimizing policy. We can define the error involved in this policy as follows.

(10) $${\rm{Error}}\left( {{\rm{minimizing}}} \right){ = ^{{\rm{df}}}} - \raise.5ex\hbox{$\scriptstyle 1$}\kern-.1em/ \kern-.15em\lower.25ex\hbox{$\scriptstyle 2$} {\rm{AIC}}\left( {{M_{min}}\left( {\bf{y}} \right)} \right) - {\rm{A}}\left( {{M_{min}}\left( {\bf{y}} \right)} \right),$$

where M _min (y) is the model with the minimum AIC score given the data y. The minimizing policy seems quite unproblematic, but here is an interesting fact: the expectation of Error(minimizing) is positive; that is, the AIC score of M _min (y) is a biased estimator of its predictive accuracy. Generally, if $\hat r$ and $\hat s$ are unbiased estimators of r and s, max( $\hat r$ , $\hat s$ ) is usually a biased estimator (in fact over-estimator) of max(r,s). Suppose r is in fact bigger than s. Obviously, max( $\hat r$ , $\hat s$ ) ≥ $\hat r$ , and after taking expectation, E(max( $\hat r$ , $\hat s$ )) ≥ E( $\hat r$ ) = r = max(r,s). However, if P( $\hat r$ < $\hat s$ ) > 0, max( $\hat r$ , $\hat s$ ) > $\hat r$ with positive probability and therefore, E(max( $\hat r$ , $\hat s$ )) > max(r, s). But unless r and s are massively different or $\hat r$ and $\hat s$ are estimators with extremely small variances, P( $\hat r$ < $\hat s$ ) > 0. In fact, in the case of AIC scores, even if model A is significantly predictively more accurate than model B, there is still a positive (however small) probability that one acquires a data set for which AIC(B) < AIC(A).

The expected value of Error(minimizing) is normally small. However, its size grows with an increase in the number of models. Without going into the technical details, we will only gesture towards an explanation for why this happens.Footnote ²⁴ Recall that the AIC score of a model, unlike its predictive accuracy, is a random variable; it varies with different data sets. Equation (4) states that the average value of AIC equals predictive accuracy, but it doesn’t say anything else about the distribution of AIC. In the most straightforward cases, −½AIC has an approximately chi-squared distribution plus a constant. Thus, in an inference problem with N models, M ₁ , M ₂ , …, M _N , we can think of −½AICs as N approximately chi-squared distributed random variables (plus a constant) with means equal to predictive accuracies, A(M ₁ ), A(M ₂ ), …, A(M _N ). For simplicity, suppose that predictive accuracies are not very different. (Neither this assumption nor the chi-square distribution is necessary for the validity of our argument, but they make things easier for our explanatory purposes.) It is unlikely for each such variable to be significantly smaller than the mean value, but if N is large, it is very likely that at least one of the variables (and thereby the minimum of all of them) is significantly smaller than the mean. Therefore, the more candidate models you have, the more biased the AIC value of M _min will be.

The fact that AIC is an unbiased estimator of predictive accuracy can mislead one into thinking that one can add to one’s candidate models at will, in the hope that if any model is not plausible it will have a poor AIC score and will thus be discarded. This is a mistaken idea. Choosing the model with minimum AIC score can be a hopelessly misguided practice if one isn’t stingy about which models to consider in the first place.

8. Concluding Remarks

We discussed Forster and Sober’s solution to the sub-family problem. Although we agree with much of what they say, we disagree about a potential implication of their argument concerning the singleton family problem. We offered our own solution for the sub-family problem, which makes the difference between the two problems salient.

Although we find the sub-family problem interesting in itself, we believe a much more important lesson about simplicity-favoring considerations can be learned from our solution to the problem. Here we would like briefly to discuss what is going on beyond the technicalities. The fundamental difference between the sub-family problem and the singleton family problem is that the sub-family model is itself a random variable dependent on and tailored to the data. (Hence the difference between equations [4] and [8].) Simple models that are designed to have excellent goodness-of-fit with the extant data tend to perform poorly in predicting future data. ^{Footnote 25} Here is a nontechnical explanation for this. Consider the very idea behind the Akaikean framework. Why is it a bad idea to use the goodness-of-fit of the best fitting member of a model as an estimate of the model’s predictive accuracy? Because to do so would essentially amount to using the current data twice: both in determining the best fitting member of the model and in determining how close the model is to the data, which is measured in terms of the fit of that same best fitting member. If the data was not used to pick out the representative (best fitting) member of the model (as is the case in singleton models), the fit of the model with the extant data was a good (unbiased) estimator of its predictive accuracy. Thus, the source of the problem with this proposal is essentially the double-use of the data. ^{Footnote 26} But the more complex the model is, the more effective such double-use will be, because the data will have more power in selecting among the members of the model. That is why the bias in taking the goodness-of-fit of the model as an estimator of its predictive accuracy increases when the complexity of the model increases. The beauty of Akaike’s results is in offering a way to calculate this bias. Now, when one designs one’s model to be simple and to have an excellent degree of fit with the current data, one re-introduces that bias into one’s estimation of predictive accuracy. In such an event, AIC is no longer an unbiased estimator, because the bias introduced by the double-use of the data is not relevant only to the number of adjustable parameters (which AIC corrects for) but also to the construction of the model itself (which AIC does not correct for).

We also discussed another problematic practice, which involves considering too many candidate models. We showed that one can effectively re-introduce the same error involved in the sub-family problem by engaging in an extreme version of this practice. Although post hoc model construction and comparing too many candidate models are problematic for two different reasons, there is a unified solution for both, namely, to respect IMR.

Before ending the paper, we would like to mention an important practical consequence of IMR. Respecting IMR gives rise to a clear criterion for determining whether a given model is gerrymandered or not. In the highly artificial examples in which the sub-family problem or similar problems are usually formulated (such as in Kukla Reference Kukla1995), it is crystal clear which models are fudged or gerrymandered (e.g., the sub-family model). However, in more realistic cases, it is sometimes not so clear. Thus, Douglas and Magnus write: “it would be perverse to do this arbitrarily, but in the general case of n-parameter models it may be possible to motivate specific values for some of the parameters. There is no formal rule for when this is or is not legitimate” (Douglas and Magnus Reference Douglas and Magnus2013, 583). A comparison between the sub-family problem and the singleton family problem suggests exactly such a rule: models in which certain parameters are held fixed are not fudged just in case there are grounds independently of the current data for holding them so fixed.