<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/checking-residual-distributions-for-non-normal-glms/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.seascapemodels.org/rstats/2017/10/06/qqplot-non-normal-glm.html"> Bluecology blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h1>Checking residual distributions for non-normal GLMs</h1>
<h2>Quantile-quantile plots</h2>
<p>If you are fitting a linear regression with Gaussian (normally<br>
distributed) errors, then one of the standard checks is to make sure the<br>
residuals are approximately normally distributed.</p>
<p>It is a good idea to do these checks for non-normal GLMs too, to make<br>
sure your residuals approximate the model’s assumption.</p>
<p>Here I explain how to create quantile-quantile plots for non-normal<br>
data, using an example of fitting a GLM using Student-t distributed<br>
errors. Such models can be appropriate when the residuals <a rel="nofollow" href="http://www.seascapemodels.org/rstats/2017/06/25/what-is-underdispersed.html" target="_blank">are<br>
overdispersed</a>.</p>
<p>First let’s create some data. We will make a linear predictor (ie the<br>
true regression line) <code>eta</code> and then simulate some data by adding<br>
residuals. We will simulate two data-sets that have the same linear<br>
predictor, but the first will have normally distributed errors and the<br>
second will have t distributed errors:</p>
<div>
<pre><code>n &lt;- 100
phi &lt;- 0.85
mu &lt;- 0.5
set.seed(23)
x &lt;- rnorm(n)

eta &lt;- mu + phi * x
nu &lt;- 2.5
tau &lt;- 3

y_tdist &lt;- eta + (rt(n, df=nu)/sqrt(tau))
y_normdist &lt;- eta + rnorm(n, sd = 1/sqrt(tau))

plot(x, y_tdist)
points(x, y_normdist, col = &quot;red&quot;, pch = 16, cex = 0.8)
legend(&quot;topleft&quot;, legend = c(&quot;t-distributed errors&quot;, &quot;normally distributed errors&quot;), pch = c(1,16), col = c(&quot;black&quot;, &quot;red&quot;))
</code></pre>
</div>
<p><img alt="" src="https://i0.wp.com/www.seascapemodels.org/images/qqplot-non-normal-glm_files/figure-markdown_strict/unnamed-chunk-1-1.png?w=456" data-recalc-dims="1"></p>
<p>Notice how the t-distributed data are more spread out. The <code>df</code><br>
parameter, here named <code>nu=2.5</code> controls how dispersed the data are.<br>
Lower values will give data that are more dispersed, large values<br>
approach a normal distribution.</p>
<p>Now let’s fit a Gaussian glm (just a linear regression really) to both<br>
these data-sets</p>
<div>
<pre><code>m1_norm &lt;- glm(y_normdist ~ x)
m1_t &lt;-  glm(y_tdist ~ x)
</code></pre>
</div>
<p>We should check whether the two models meet the normal assumption, using<br>
the standard ‘qq’ (quantile-quantile) plot:</p>
<div>
<pre><code>par(mfrow = c(1,2))
plot(m1_norm, 2)
plot(m1_t, 2)
</code></pre>
</div>
<p><img alt="" src="https://i2.wp.com/www.seascapemodels.org/images/qqplot-non-normal-glm_files/figure-markdown_strict/unnamed-chunk-3-1.png?w=456" data-recalc-dims="1"></p>
<p>These plots compare the theoretical quantiles to the actual quantiles of<br>
the residuals. If the points fall on the straight line then the<br>
theoretical and realised are very similar, and the assumption is met.<br>
Clearly this is not the case for the second model, which is<br>
overdispersed.</p>
<p>We know it is overdispersed because the theoretical quantiles are much<br>
smaller than the actual at the tails (notice how the ends down then up).</p>
<p>The p-values (or CIs if you use them) for <code>m1_t</code> are therefore likely<br>
biased and too narrow, leading potentially to type I errors (us saying<br>
that x affects y, which in fact it does not). In this case we know we<br>
haven’t made a type I error, because we made up the data. However, if<br>
you were using real data you wouldn’t be so sure.</p>
<h2>Doing our own quantile-quantile plot</h2>
<p>To better understand the QQ plot it helps to generate it yourself,<br>
rather than using R’s automatic checks.</p>
<p>First we calculate the model residuals (in <code>plot(m1_t)</code> R did this<br>
internally):</p>
<div>
<pre><code>m1_t_resid &lt;- y_tdist - predict(m1_t)
</code></pre>
</div>
<p>Then we can plot the quantiles for the residuals against theoretical<br>
quantiles generated using <code>qnorm</code>. Below we also plot the original QQ<br>
plot from above, so you can see that our version is the same as R’s<br>
automatic one:</p>
<div>
<pre><code>par(mfrow = c(1,2))
qqplot(qnorm(ppoints(n), sd = 1), m1_t_resid)
qqline(m1_t_resid, lty = 2)

plot(m1_t,2)
</code></pre>
</div>
<p><img alt="" src="https://i2.wp.com/www.seascapemodels.org/images/qqplot-non-normal-glm_files/figure-markdown_strict/unnamed-chunk-5-1.png?w=456" data-recalc-dims="1"></p>
<p>I added the <code>qqline</code> for comparative purposes. It just puts a line<br>
through the 25th and 75th quantiles.</p>
<h2>QQ plot for a non-normal GLM</h2>
<p>Now we have learned how to write our own custom for a QQ plot, we can<br>
use it to check other types of non-normal data.</p>
<p>Here we will fit a GLM to the <code>y_tdist data</code> using student-t distributed<br>
errors. I do this using the Bayesian package <code>INLA</code>.</p>
<div>
<pre><code>library(INLA)

data &lt;- list(y=y_tdist, x = x)
mod_tdist &lt;- inla(y ~ x, family=&quot;T&quot;, data=data,
    control.predictor = list(compute = TRUE),
    control.family =
        list(
            hyper =
                list(prec = list(prior=&quot;loggamma&quot;,param=c(1,0.5)),
                    dof = list(prior=&quot;loggamma&quot;,param=c(1,0.5))
                    )
            )
        )
</code></pre>
</div>
<p>The <code>family =&quot;T&quot;</code> command tells INLA to use the t-distribution, rather<br>
than the Normal distribution. Note also I have specified the priors<br>
using the <code>control.family</code> command. This is best practice. We need a<br>
prior for the precision (1/variance) and a prior for the <code>dof</code> (=<br>
degrees of freedom, which has to be >2 in INLA).</p>
<p>It is sometimes help to visualise the priors, so we can check too see<br>
they look sensible. Here we visualise the prior for the <code>dof</code>, (which in<br>
INLA has a min value of 2):</p>
<div>
<pre><code>xgamma &lt;- seq(0.01, 10, length.out = 100)
plot(xgamma+2, dgamma(xgamma, 1, 0.5), type = 'l', xlab = &quot;Quantile&quot;, ylab = &quot;Density&quot;)
</code></pre>
</div>
<p><img alt="" src="https://i0.wp.com/www.seascapemodels.org/images/qqplot-non-normal-glm_files/figure-markdown_strict/unnamed-chunk-7-1.png?w=456" data-recalc-dims="1"></p>
<p>We don’t really expect values much greater than 10, so this prior makes<br>
sense. If we used an old-school prior that was flat in 2-1000 we might<br>
get issues with model fitting.</p>
<p>Now enough about priors. Let’s look at the estimated coefficients:</p>
<div>
<pre><code>mod_tdist$summary.fixed

##                  mean         sd 0.025quant  0.5quant 0.975quant      mode
## (Intercept) 0.5324814 0.07927198  0.3773399 0.5321649  0.6891779 0.5315490
## x           0.7229362 0.08301006  0.5565746 0.7239544  0.8835630 0.7259817
##                      kld
## (Intercept) 3.067485e-12
## x           6.557565e-12
</code></pre>
</div>
<p>Good the CIs contain our true values, and the mean is close to our true<br>
value too.<br>
What about the hyper-parameters (the precision and DF)? We need to get<br>
INLA to run some more calucations to get accurate estimates of these:</p>
<div>
<pre><code>h_tdist &lt;- inla.hyperpar(mod_tdist)
h_tdist$summary.hyperpar[,3:5]

##                                          0.025quant  0.5quant 0.975quant
## precision for the student-t observations  0.2663364 0.6293265   1.163440
## degrees of freedom for student-t          2.2404966 2.7396391   4.459057
</code></pre>
</div>
<p>The estimate for the DF might be a ways off the mark. That is ok, we<br>
expect that, you need lots of really good data to get accurate estimates<br>
of hyper-parameters.</p>
<p>Now, let’s use our skills in creating QQ plots to make QQ plot using<br>
theoretical quantiles from the t distribution.</p>
<p>First step is to extract INLA’s predictions of the data, so we can<br>
calculate residuals</p>
<div>
<pre><code>preds &lt;- mod_tdist$summary.fitted.values
resids &lt;- y_tdist - preds[,4]
</code></pre>
</div>
<p>Next step is to extract the marginal estimates of the DF and precision<br>
to use when generating our QQ plot (the quantiles will change with the<br>
DF):</p>
<div>
<pre><code>tau_est &lt;- h_tdist$summary.hyperpar[1,4]
nu_est &lt;- h_tdist$summary.hyperpar[2,4]
</code></pre>
</div>
<p>Now we can use <code>qt()</code> to generate theoretical quantiles and the<br>
residuals for our realised quantiles:</p>
<div>
<pre><code>qqplot(qt(ppoints(n), df = nu_est), resids * sqrt(tau_est),
    xlab = &quot;Theoretical quantile&quot;, ylab = &quot;residuals&quot;)
qqline(resids * sqrt(tau_est), lty = 2)
</code></pre>
</div>
<p><img alt="" src="https://i0.wp.com/www.seascapemodels.org/images/qqplot-non-normal-glm_files/figure-markdown_strict/unnamed-chunk-12-1.png?w=456" data-recalc-dims="1"></p>
<p>Note that I multiply the residuals by the sqrt of the precision<br>
estimate. This is <a rel="nofollow" href="http://www.math.ntnu.no/inla/r-inla.org/doc/likelihood/student-t.pdf" target="_blank">how INLA fits a t-distributed<br>
GLM</a>.<br>
I do the same for the <code>qqline</code>.</p>
<p>Our residuals are now falling much closer to the line. The model is<br>
doing a much better job of fitting the data. You could also calculate<br>
the WAIC for this model and a Gaussian one, to compare the fits. The<br>
t-distributed GLM should have a lower WAIC (better fit).</p>
<p>We can now be confident that our CIs are accurate.</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/2017/10/06/qqplot-non-normal-glm.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/fj0boprtXH8" height="1"></a></div>