Repository: statsmodels/statsmodels (Python)

Issue:
AIC value inconsistent with R
#### Describe the bug

The definition of AIC in linear regression (sm.api.OLS()) does not match that of R. This is due to the fact that the number of parameters is defined as (df_model + 1), which does not take into account the variance term that is also estimated from the data. (This is why this bug happens in OLS but not in case of GLM.) I think this should have been (df_model + 2). This bug actually does not change rankings of models because the +1 term vanishes when comparing models, but the definition is conceptually wrong and would potentially cause incompatibilities with some of R functions.

#### Code Sample, a copy-pastable example if possible

import statsmodels.api as sm
from sklearn.datasets import load_iris
from numpy import c_, ones, log

data = load_iris()

X = c_[ones(data['data'].shape[0]), data['data'][:, 1]]
y = data['data'][:, 0]
n = data['data'].shape[0]

model = sm.OLS(y, X)
fit = model.fit()

print(fit.aic) # gives 369.9916713254629
print(fit.bic) # gives 376.0129419136554

print(-2*fit.llf + 3*2) # R AIC: 371.9916713254629
print(-2*fit.llf + 3*log(n)) # R BIC: 381.02357720775166

```python
# Your code here that produces the bug
# This example should be self-contained, and so not rely on external data.
# It should run in a fresh ipython session, and so include all relevant imports.
```
<details>

**Note**: As you can see, there are many issues on our GitHub tracker, so it is very possible that your issue has been posted before. Please check first before submitting so that we do not have to handle and close duplicates.

**Note**: Please be sure you are using the latest released version of `statsmodels`, or a recent build of `master`. If your problem has been fixed in an unreleased version, you might be able to use `master` until a new release occurs. 

**Note**: If you are using a released version, have you verified that the bug exists in the master branch of this repository? It helps the limited resources if we know problems exist in the current master so that they do not need to check whether the code sample produces a bug in the next release.

</details>


If the issue has not been resolved, please file it in the issue tracker.

#### Expected Output

A clear and concise description of what you expected to happen.

#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

[paste the output of ``import statsmodels.api as sm; sm.show_versions()`` here below this line]

</details>



Hints:
same issue for BIC #1802

no decision to change it yet
Label this as a FAQ. I think we are past the point where this can change.
Still showing in 0.12
If it's wrong, could it be corrected?
It's not just a diagnostic: if you're using the AIC/BIC to compare models, it can have a material effect.

It does _not_ affect the relative rank of models with the same IC which is how IC are used in model selection.  
The plan is to add an `infocrit` method to results classes that has options. GLM also does not count an estimated scale or dispersion parameter in the ic.

In OLS (and GLM), the parameter estimates and fit in terms of sum of squares, or residual measures is independent of a scale estimate. So there is not necessarily a reason to include the scale estimate when comparing models with different exog.

But we do want consistent definitions across models, so we can us the ic for choosing among those, e.g. count models with different distributional assumptions.



Just a note that state space models have an `info_criteria` method rather than an `infocrit` method.
@ChadFulton Thanks, then we can coordinate on that.
I haven't looked at tsa in some time.)
@ChadFulton 
Is scale included in the count of parameters for ic in tsa and statespace?
> Is scale included in the count of parameters for ic in tsa and statespace?

TL;DR: `SARIMAX`, `ARIMA`, and `AutoReg` all include the scale in the definition of AIC, while `VAR` does not.

--------

Scale is included in `df_model` and in the information criteria computations for state space models (e.g. `SARIMAX`). This is true both when it is estimated as part of the overall parameter vector and when we use the version where we first concentrate it out of the likelihood.

I took a look at the old `ARIMA` model, and there it is also included in the information criteria. This model always estimates the scale by concentrating it out, however, `df_model` does not include report it (so the information criteria computations use the private `_ic_df_model`.

The IC for both `AutoReg` and `VAR` are based on Lutkepohl's definition, and they (a) do not include the intercept terms, (b) drops the constant term, and (c) divide through by `nobs` (in fact, [Stata's `varsoc`](https://www.stata.com/manuals13/tsvarsoc.pdf) divides even its standard AIC through by `nobs`). However, `AutoReg` does include the scale term as a parameter (different from Lutkepohl's formulation) while `VAR` does not include the free terms in the error covariance matrix (the same as Lutkepohl's formulation).

--------

One other thing that is interesting / annoying: the way that Lutkepohl forms the AIC by concentrating out the error covariance matrix implies that the constant term that is dropped is of the form `k_endog * (np.log(2 * np.pi) + 1)`, while in other cases (such as R's `forecast` package, if I remember correctly), the constant term that is dropped is of the form `k_endog * np.log(2 * np.pi)`. (Where, again, either of these cases may or may not also be multipled by `nobs`, depending on the definition used).
@ChadFulton Thanks for the overview

sounds also pretty messy outside full MLE

One problem is that in the least squares and quasi-likelihood version we don't need a fully specified likelihood and variance/scale estimate are in most of those cases orthogonal to the mean parameter estimates.

The traditional sum of squares, Sigma definition of aic divides by nobs
https://www.statsmodels.org/devel/generated/statsmodels.tools.eval_measures.aic_sigma.html

I'm only part way to understanding more of the IC for not fully specified or misspecified likelihood cases.
e.g. TIC and GIC #6516
I still need to see how this works out in various models with specified or un-/misspecified likelihood function.

Adding options for various definitions of IC pushes the decision which to use to the user :)

As Kevin pointed out, any constant term is irrelevant for specification decisions within a model. But if we have different models or different distributional assumptions for the same underlying modelling task, then we need consistent definitions, whatever they are.

Another issue: What is relevant in the IC might also depend on the target of the estimation.
In OLS and similar we are interested in the mean function and we might not care about the exact distribution of the residuals. So scale and higher order properties of OLS or Poisson model are not of interest.
If we are interested in a correctly specified distribution, then we don't want to ignore those extra features of the distribution.
#7142

(I guess a tsa example would be an ARMA-GARCH that is separable in mean and variance model. For the specification search for the mean model, we don't really need to include all parameters from the GARCH model. If we want to check the overall fit, then we need the full model.
My similar examples these days are count models and overdispersion)


Patch:
diff --git a/statsmodels/genmod/generalized_linear_model.py b/statsmodels/genmod/generalized_linear_model.py
--- a/statsmodels/genmod/generalized_linear_model.py
+++ b/statsmodels/genmod/generalized_linear_model.py
@@ -1768,7 +1768,7 @@ def aic(self):
         Akaike Information Criterion
         -2 * `llf` + 2 * (`df_model` + 1)
         """
-        return -2 * self.llf + 2 * (self.df_model + 1)
+        return self.info_criteria("aic")
 
     @property
     def bic(self):
@@ -1826,11 +1826,9 @@ def bic_llf(self):
         Based on the log-likelihood,
         -2 * `llf` + log(n) * (`df_model` + 1)
         """
-        return -2*self.llf + (self.df_model+1)*np.log(
-            self.df_model+self.df_resid+1
-        )
+        return self.info_criteria("bic")
 
-    def info_criteria(self, crit, scale=None):
+    def info_criteria(self, crit, scale=None, dk_params=0):
         """Return an information criterion for the model.
 
         Parameters
@@ -1840,6 +1838,11 @@ def info_criteria(self, crit, scale=None):
         scale : float
             The scale parameter estimated using the parent model,
             used only for qaic.
+        dk_params : int or float
+            Correction to the number of parameters used in the information
+            criterion. By default, only mean parameters are included, the
+            scale parameter is not included in the parameter count.
+            Use ``dk_params=1`` to include scale in the parameter count.
 
         Returns the given information criterion value.
 
@@ -1860,13 +1863,15 @@ def info_criteria(self, crit, scale=None):
         Burnham KP, Anderson KR (2002). Model Selection and Multimodel
         Inference; Springer New York.
         """
-
         crit = crit.lower()
+        k_params = self.df_model + 1 + dk_params
 
         if crit == "aic":
-            return self.aic
+            return -2 * self.llf + 2 * k_params
         elif crit == "bic":
-            return self.bic
+            nobs = self.df_model + self.df_resid + 1
+            bic = -2*self.llf + k_params*np.log(nobs)
+            return bic
         elif crit == "qaic":
             f = self.model.family
             fl = (families.Poisson, families.NegativeBinomial,
@@ -1876,7 +1881,7 @@ def info_criteria(self, crit, scale=None):
                 msg += "Negative Binomial families."
                 warnings.warn(msg)
             llf = self.llf_scaled(scale=1)
-            return -2 * llf/scale + 2 * (self.df_model + 1)
+            return -2 * llf/scale + 2 * k_params
 
     @Appender(pred.get_prediction_glm.__doc__)
     def get_prediction(self, exog=None, exposure=None, offset=None,
diff --git a/statsmodels/regression/linear_model.py b/statsmodels/regression/linear_model.py
--- a/statsmodels/regression/linear_model.py
+++ b/statsmodels/regression/linear_model.py
@@ -1826,7 +1826,7 @@ def aic(self):
         For a model with a constant :math:`-2llf + 2(df\_model + 1)`. For a
         model without a constant :math:`-2llf + 2(df\_model)`.
         """
-        return -2 * self.llf + 2 * (self.df_model + self.k_constant)
+        return self.info_criteria("aic")
 
     @cache_readonly
     def bic(self):
@@ -1836,8 +1836,42 @@ def bic(self):
         For a model with a constant :math:`-2llf + \log(n)(df\_model+1)`.
         For a model without a constant :math:`-2llf + \log(n)(df\_model)`.
         """
-        return (-2 * self.llf + np.log(self.nobs) * (self.df_model +
-                                                     self.k_constant))
+        return self.info_criteria("bic")
+
+    def info_criteria(self, crit, dk_params=0):
+        """Return an information criterion for the model.
+
+        Parameters
+        ----------
+        crit : string
+            One of 'aic', 'bic', 'aicc' or 'hqic'.
+        dk_params : int or float
+            Correction to the number of parameters used in the information
+            criterion. By default, only mean parameters are included, the
+            scale parameter is not included in the parameter count.
+            Use ``dk_params=1`` to include scale in the parameter count.
+
+        Returns the given information criterion value.
+
+        References
+        ----------
+        Burnham KP, Anderson KR (2002). Model Selection and Multimodel
+        Inference; Springer New York.
+        """
+        crit = crit.lower()
+        k_params = self.df_model + self.k_constant + dk_params
+
+        if crit == "aic":
+            return -2 * self.llf + 2 * k_params
+        elif crit == "bic":
+            bic = -2*self.llf + np.log(self.nobs) * k_params
+            return bic
+        elif crit == "aicc":
+            from statsmodels.tools.eval_measures import aicc
+            return aicc(self.llf, self.nobs, k_params)
+        elif crit == "hqic":
+            from statsmodels.tools.eval_measures import hqic
+            return hqic(self.llf, self.nobs, k_params)
 
     @cache_readonly
     def eigenvals(self):
diff --git a/statsmodels/regression/rolling.py b/statsmodels/regression/rolling.py
--- a/statsmodels/regression/rolling.py
+++ b/statsmodels/regression/rolling.py
@@ -525,6 +525,10 @@ def bic(self):
         with np.errstate(divide="ignore"):
             return self._wrap(RegressionResults.bic.func(self))
 
+    def info_criteria(self, crit, dk_params=0):
+        return self._wrap(RegressionResults.info_criteria(
+            self, crit, dk_params=dk_params))
+
     @cache_readonly
     def params(self):
         """Estimated model parameters"""


Test patch:
diff --git a/statsmodels/genmod/tests/test_glm.py b/statsmodels/genmod/tests/test_glm.py
--- a/statsmodels/genmod/tests/test_glm.py
+++ b/statsmodels/genmod/tests/test_glm.py
@@ -2540,6 +2540,7 @@ def test_output_exposure_null(reset_randomstate):
     # Check that they are different
     assert np.abs(null_model_without_exposure.llf - model.llnull) > 1
 
+
 def test_qaic():
 
     # Example from documentation of R package MuMIn
@@ -2562,7 +2563,11 @@ def test_qaic():
     # presumably because they count the scale parameter in df.
     # This won't matter when comparing models by differencing
     # QAICs.
+    # Binomial doesn't have a scale parameter, so adding +1 is not correct.
     assert_allclose(qaic, 29.13266, rtol=1e-5, atol=1e-5)
+    qaic1 = r.info_criteria(crit="qaic", scale=scale, dk_params=1)
+    assert_allclose(qaic1, 31.13266, rtol=1e-5, atol=1e-5)
+
 
 def test_tweedie_score():
 
diff --git a/statsmodels/regression/tests/test_regression.py b/statsmodels/regression/tests/test_regression.py
--- a/statsmodels/regression/tests/test_regression.py
+++ b/statsmodels/regression/tests/test_regression.py
@@ -173,6 +173,15 @@ def test_loglike(self):
 
     def test_aic(self):
         assert_almost_equal(self.res1.aic, self.res2.aic, self.decimal_aic)
+        # the following just checks the definition
+        aicc1 = self.res1.info_criteria("aicc")
+        k = self.res1.df_model + self.res1.model.k_constant
+        nobs = self.res1.model.nobs
+        aicc2 = self.res1.aic + 2 * (k**2 + k) / (nobs - k - 1)
+        assert_allclose(aicc1, aicc2, rtol=1e-10)
+        hqic1 = self.res1.info_criteria("hqic")
+        hqic2 = (self.res1.aic - 2 * k) + 2 * np.log(np.log(nobs)) * k
+        assert_allclose(hqic1, hqic2, rtol=1e-10)
 
     decimal_bic = DECIMAL_4