<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/how-do-i-interpret-the-aic/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.seascapemodels.org/rstats/2018/04/13/how-to-use-the-AIC.html"> Bluecology blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h1>How do I interpret the AIC?</h1>
<p>My student asked today how to interpret the AIC (Akaike’s Information<br>
Criteria) statistic for model selection. We ended up bashing out some R<br>
code to demonstrate how to calculate the AIC for a simple GLM (general<br>
linear model). I always think if you can understand the derivation of a<br>
statistic, it is much easier to remember how to use it.</p>
<p>Now if you google derivation of the AIC, you are likely to run into a<br>
lot of math. <strong>But the principles are really not that complex.</strong> So here<br>
we will fit some simple GLMs, then derive a means to choose the ‘best’<br>
one.</p>
<p>Skip to the end if you just want to go over the basic principles.</p>
<p>Before we can understand the AIC though, we need to understand the<br>
statistical methodology of likelihoods.</p>
<h2>Explaining likelihoods</h2>
<p>Say you have some data that are normally distributed with a mean of <code>5</code><br>
and an sd of <code>3</code>:</p>
<div>
<div>
<pre><code>set.seed(126)
n &lt;- 50 #sample size
a &lt;- 5
sdy &lt;- 3
y &lt;- rnorm(n, mean = a, sd = sdy)
hist(y)
</code></pre>
</div>
</div>
<p><img alt="" src="https://i0.wp.com/www.seascapemodels.org/images/how-to-use-the-AIC_files/figure-markdown_strict/unnamed-chunk-1-1.png?w=456" data-recalc-dims="1"></p>
<p>Now we want to estimate some parameters for the population that <code>y</code> was<br>
sampled from, like its mean and standard devaiation (which we know here<br>
to be 5 and 3, but in the real world you won’t know that).</p>
<p>We are going to use frequentist statistics to estimate those parameters.<br>
Philosophically this means we believe that there is ‘one true value’ for<br>
each parameter, and the data we observed are generated by this true<br>
value.</p>
<div>
<div>
<pre><code>m1 &lt;- glm(y ~ 1, family = &quot;gaussian&quot;)
sm1 &lt;- summary(m1)
</code></pre>
</div>
</div>
<p>The estimate of the mean is stored here <code>coef(m1)</code> =4.38, the estimated<br>
variance here <code>sm1$dispersion</code>= 5.91, or the SD <code>sqrt(sm1$dispersion)</code><br>
=2.43. Just to be totally clear, we also specified that we believe the<br>
data follow a normal (AKA “Gaussian”) distribution.</p>
<p>We just fit a GLM asking R to estimate an intercept parameter (<code>~1</code>),<br>
which is simply the mean of <code>y</code>. We also get out an estimate of the SD<br>
(= $\sqrt variance$) You might think its overkill to use a GLM to<br>
estimate the mean and SD, when we could just calculate them directly.</p>
<p>Well notice now that R also estimated some other quantities, like the<br>
residual deviance and the AIC statistic.</p>
<div>
<div>
<pre><code>summary(m1)

##
## Call:
## glm(formula = y ~ 1, family = &quot;gaussian&quot;)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -5.7557  -0.9795   0.2853   1.7288   3.9583
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   4.3837     0.3438   12.75   &lt;2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 5.910122)
##
##     Null deviance: 289.6  on 49  degrees of freedom
## Residual deviance: 289.6  on 49  degrees of freedom
## AIC: 233.72
##
## Number of Fisher Scoring iterations: 2
</code></pre>
</div>
</div>
<p>You might also be aware that the deviance is a measure of model fit,<br>
much like the sums-of-squares. Note also that the value of the AIC is<br>
suspiciously close to the deviance. Despite its odd name, the concepts<br>
underlying the deviance are quite simple.</p>
<p>As I said above, we are observing data that are generated from a<br>
population with one true mean and one true SD. Given we know have<br>
estimates of these quantities that define a probability distribution, we<br>
could also estimate the likelihood of measuring a new value of <code>y</code> that<br>
say = 7.</p>
<p>To do this, we simply plug the estimated values into the equation for<br>
the normal distribution and ask for the relative likelihood of <code>7</code>. We<br>
do this with the R function <code>dnorm</code></p>
<div>
<div>
<pre><code>sdest &lt;- sqrt(sm1$dispersion)
dnorm(7, mean = coef(m1), sd = sdest)

## [1] 0.09196167
</code></pre>
</div>
</div>
<p>Formally, this is the relative likelihood of the value 7 given the<br>
values of the mean and the SD that we estimated (=4.8 and 2.39<br>
respectively if you are using the same random seed as me).</p>
<p>You might ask why the likelihood is greater than 1, surely, as it comes<br>
from a probability distribution, it should be &lt;1. Well, the normal<br>
distribution is continuous, which means it describes an infinte set of<br>
possible <code>y</code> values, so the probability of any given value will be zero.<br>
The relative likelihood on the other hand can be used to <a rel="nofollow" href="https://en.wikipedia.org/wiki/Probability_density_function" target="_blank">calculate the<br>
probability of a range of<br>
values</a>.</p>
<p>So you might realise that calculating the likelihood of all the data<br>
would be a sensible way to measure how well our ‘model’ (just a mean and<br>
SD here) fits the data.</p>
<p>Here’s what the likelihood looks like:</p>
<div>
<div>
<pre><code>plot(y, dnorm(y, mean = coef(m1), sd = sdest), ylab = &quot;Likelihood&quot;)
</code></pre>
</div>
</div>
<p><img alt="" src="https://i0.wp.com/www.seascapemodels.org/images/how-to-use-the-AIC_files/figure-markdown_strict/unnamed-chunk-5-1.png?w=456" data-recalc-dims="1"></p>
<p>It’s just a normal distribution.</p>
<p>To do this, think about how you would calculate the probability of<br>
multiple (independent) events. Say the chance I ride my bike to work on<br>
any given day is 3/5 and the chance it rains is 161/365 (like<br>
Vancouver!), then the chance I will ride in the rain[1] is 3/5 *<br>
161/365 = about 1/4, so I best wear a coat if riding in Vancouver.</p>
<p>We can do the same for likelihoods, simply multiply the likelihood of<br>
each individual <code>y</code> value and we have the total likelihood. This will be<br>
a very small number, because we multiply a lot of small numbers by each<br>
other. So one trick we use is to sum the log of the likelihoods instead<br>
of multiplying them:</p>
<div>
<div>
<pre><code>y_lik &lt;- dnorm(y, mean = coef(m1), sd = sdest, log = TRUE)
sum(y_lik)

## [1] -114.8636
</code></pre>
</div>
</div>
<p>The larger (the less negative) the likelihood of our data given the<br>
model’s estimates, the ‘better’ the model fits the data. The deviance is<br>
calculated from the likelihood and for the deviance smaller values<br>
indicate a closer fit of the model to the data.</p>
<p>The parameter values that give us the smallest value of the<br>
-log-likelihood are termed the maximum likelihood estimates.</p>
<h2>Comparing alternate hypotheses with likelihoods</h2>
<p>Now say we have measurements and two covariates, <code>x1</code> and <code>x2</code>, either<br>
of which we think might affect y:</p>
<div>
<div>
<pre><code>a &lt;- 5
b &lt;- 3
n &lt;- 100
x1 &lt;- rnorm(n)
x2 &lt;- rnorm(n)
sdy &lt;- 1
y &lt;- a + b*x1 + rnorm(n, sd = sdy)
par(mfrow = c(1,2))
plot(x1, y)
plot(x2, y)
</code></pre>
</div>
</div>
<p><img alt="" src="https://i2.wp.com/www.seascapemodels.org/images/how-to-use-the-AIC_files/figure-markdown_strict/unnamed-chunk-7-1.png?w=456" data-recalc-dims="1"></p>
<p>So x1 is a cause of y, but x2 does not affect y. How would we choose<br>
which hypothesis is most likely? Well one way would be to compare models<br>
with different combinations of covariates:</p>
<div>
<div>
<pre><code>m1 &lt;- glm(y ~ x1)
m2 &lt;- glm(y ~ x2)
m3 &lt;- glm(y ~ x1 + x2)
</code></pre>
</div>
</div>
<p>Now we are fitting a line to y, so our estimate of the mean is now the<br>
line of best fit, it varies with the value of x1. To visualise this:</p>
<div>
<div>
<pre><code>plot(x1, y)
lines(x1, predict(m1))
</code></pre>
</div>
</div>
<p><img alt="" src="https://i0.wp.com/www.seascapemodels.org/images/how-to-use-the-AIC_files/figure-markdown_strict/unnamed-chunk-9-1.png?w=456" data-recalc-dims="1"></p>
<p>The <code>predict(m1)</code> gives the line of best fit, ie the mean value of y<br>
given each x1 value. We then use predict to get the likelihoods for each<br>
model:</p>
<div>
<div>
<pre><code>sm1 &lt;- summary(m1)
sum(dnorm(y, mean = predict(m1), sd = sqrt(sm1$dispersion), log = TRUE))

## [1] -125.6214

sm2 &lt;- summary(m2)
sum(dnorm(y, mean = predict(m2), sd = sqrt(sm2$dispersion), log = TRUE))

## [1] -247.8059

sm3 &lt;- summary(m3)
sum(dnorm(y, mean = predict(m3), sd = sqrt(sm3$dispersion), log = TRUE))

## [1] -125.4843
</code></pre>
</div>
</div>
<p>The likelihood of <code>m1</code> is larger than <code>m2</code>, which makes sense because<br>
<code>m2</code> has the ‘fake’ covariate in it. The likelihood for <code>m3</code> (which has<br>
both x1 and x2 in it) is fractionally larger than the likelihood <code>m1</code>,<br>
so should we judge that model as giving nearly as good a representation<br>
of the data?</p>
<p>Because the likelihood is only a tiny bit larger, the addition of <code>x2</code><br>
has only explained a tiny amount of the variance in the data. But where<br>
do you draw the line between including and excluding x2? You run into a<br>
similar problem if you use R^2 for model selection.</p>
<p>So what if we penalize the likelihood by the number of paramaters we<br>
have to estimate to fit the model? Then if we include more covariates<br>
(and we estimate more slope parameters) only those that account for a<br>
lot of the variation will overcome the penalty.</p>
<p>What we want a statistic that helps us select the most parsimonious<br>
model.</p>
<h2>The AIC as a measure of parsimony</h2>
<p>One way we could penalize the likelihood by the number of parameters is<br>
to add an amount to it that is proportional to the number of parameters.<br>
First, let’s multiply the log-likelihood by -2, so that it is positive<br>
and smaller values indicate a closer fit.</p>
<div>
<div>
<pre><code>LLm1 &lt;- sum(dnorm(y, mean = predict(m1), sd = sqrt(sm1$dispersion), log = TRUE))
-2*LLm1

## [1] 251.2428
</code></pre>
</div>
</div>
<p>Why its -2 not -1, I can’t quite remember, but I think just <a rel="nofollow" href="https://en.wikipedia.org/wiki/Akaike_information_criterion#CITEREFAkaike1985" target="_blank">historical<br>
reasons</a>.</p>
<p>Then add 2*k, where k is the number of estimated parameters.</p>
<div>
<div>
<pre><code>-2*LLm1 + 2*3

## [1] 257.2428
</code></pre>
</div>
</div>
<p>For <code>m1</code> there are three parameters, one intercept, one slope and one<br>
standard deviation. Now, let’s calculate the AIC for all three models:</p>
<div>
<div>
<pre><code>-2*LLm1 + 2*3

## [1] 257.2428

LLm2 &lt;- sum(dnorm(y, mean = predict(m2), sd = sqrt(sm2$dispersion), log = TRUE))
-2*LLm2 + 2*3

## [1] 501.6118

LLm3 &lt;- sum(dnorm(y, mean = predict(m3), sd = sqrt(sm3$dispersion), log = TRUE))
-2*LLm3 + 2*4

## [1] 258.9686
</code></pre>
</div>
</div>
<p>We see that model 1 has the lowest AIC and therefore has the most<br>
parsimonious fit. Model 1 now outperforms model 3 which had a slightly<br>
higher likelihood, but because of the extra covariate has a higher<br>
penalty too.</p>
<h2>AIC basic principles</h2>
<p>So to summarize, the basic principles that guide the use of the AIC are:</p>
<ol>
<li>
<p>Lower indicates a more parsimonious model, relative to a model fit<br>
with a higher AIC.</p>
</li>
<li>
<p>It is a <em>relative</em> measure of model parsimony, so it only has<br>
meaning if we compare the AIC for alternate hypotheses (= different<br>
models of the data).</p>
</li>
<li>
<p>We can compare non-nested models. For instance, we could compare a<br>
linear to a non-linear model.</p>
</li>
<li>
<p>The comparisons are only valid for models that are fit to the same response<br>
data (ie values of y).</p>
</li>
<li>
<p>Model selection conducted with the AIC will choose the same model as<br>
leave-one-out cross validation (where we leave out one data point<br>
and fit the model, then evaluate its fit to that point) for large<br>
sample sizes.</p>
</li>
<li>
<p>You shouldn’t compare <em>too</em> many models with the AIC. You will run<br>
into the same problems with multiple model comparison as you would<br>
with p-values, in that you might by chance find a model with the<br>
lowest AIC, that isn’t truly the most appropriate model.</p>
</li>
<li>
<p>When using the AIC you might end up with multiple models that<br>
perform similarly to each other. So you have similar evidence<br>
weights for different alternate hypotheses. In the example above m3<br>
is actually about as good as m1.</p>
</li>
<li>
<p>You should correct for small sample sizes if you use the AIC with<br>
small sample sizes, by using the AICc statistic.</p>
</li>
</ol>
<p>[1] Assuming it rains all day, which is reasonable for Vancouver.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.seascapemodels.org/rstats/2018/04/13/how-to-use-the-AIC.html"> Bluecology blog</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/3Yzz6JFB2_o" height="1"></a></div>