<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/confidence-intervals-for-glms/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://www.fromthebottomoftheheap.net/2018/12/10/confidence-intervals-for-glms/"> From the Bottom of the Heap - R</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>
You’ve estimated a GLM or a related model (GLMM, GAM, etc.) for your latest paper and, like a good researcher, you want to visualise the model and show the uncertainty in it. In general this is done using confidence intervals with typically 95% converage. If you remember a little bit of theory from your stats classes, you may recall that such an interval can be produced by adding to and subtracting from the fitted values 2 times their standard error. Unfortunately this only really works like this for a linear model. If I had a dollar (even a Canadian one) for every time I’ve seen someone present graphs of estimated abundance of some species where the confidence interval includes negative abundances, I’d be rich! Here, following the rule of &quot;if I’m asked more than once I should write a blog post about it!&quot; I’m going to show a simple way to correctly compute a confidence interval for a GLM or a related model.
</p>
<h3>
Why is plus/minus two standard errors wrong?<br>
</h3>
<p>
Well, it’s not! However, the main reason why people mess up computing confidence intervals for a GLM is that they do all the calculations on the <em>response</em> scale. This results in symmetric intervals on this scale and the very real possibility that the intervals will include values that are nonsensical, like negative abundances and concentrations, or probabilities that are outside the limits of 0 and 1.
</p>
<p>
Think about a Poisson GLM fitted to some species abundance data. In this model there is an implied mean-variance relationship; as the mean count increases so does the variance. In fact, in the Poisson GLM, the mean and variance are the same thing. The implication of this is that as the mean tends to zero, so must the variance. If we had an expected count of zero the variance would also be zero, and our uncertainty about this value would also be zero. However, our model won’t ever return expected (fitted) values that are exactly equal to zero; it might yield values that are very close to zero, but never exactly zero. In that case we do have some uncertainty about this fitted value; the uncertainty on the lower end has to logically fit somewhere between the small estimated value and zero, but not exactly zero as we’re not creating an interval with 100% coverage.
</p>
<p>
We might also logically expect greater uncertainty above the fitted value, for our upper limit on the confidence interval; we’re saying that the true expected abudance is possibly somewhat larger than the fitted value and due to the mean-variance relationship, a larger fitted value is a larger mean value, which implies a larger variance, and consequently a larger amount of uncertainty above the fitted value than below.
</p>
<p>
Similar arguments can be made for models where there are both upper and lower limits to the response, such as binomial models where the response is a probability bounded between 0 and 1. As the fitted value approaches either boundary the uncertainty about the fitted value in the direction of the boundary gets squished up and the asymmetry of the confidence interval increases.
</p>
<p>
To illustrate, I’ll use a simple data set on wasp visits to leaves of the Cobra Lily, <em>Darlingtonia californica</em>. The data are on my blog and I’ve created a short link using bitly.com. If you want to follow along, load the data and some packages as shown
</p>
<figure>
<pre><code data-lang="r"><span>## packages
</span><span>library</span><span>(</span><span>'readr'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'tibble'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'dplyr'</span><span>)</span><span>
</span><span>library</span><span>(</span><span>'ggplot2'</span><span>)</span><span>
</span><span>theme_set</span><span>(</span><span>theme_bw</span><span>())</span><span>
</span><span>wasp</span><span> </span><span>&lt;-</span><span> </span><span>read_csv</span><span>(</span><span>'http://bit.ly/cobralily'</span><span>,</span><span> </span><span>skip</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span><span>wasp</span><span> </span><span>&lt;-</span><span> </span><span>mutate</span><span>(</span><span>wasp</span><span>,</span><span> </span><span>lvisited</span><span> </span><span>=</span><span> </span><span>as.logical</span><span>(</span><span>visited</span><span>))</span></code></pre>
</figure>
<p>
The experiment used timed census of visitations by wasps to leaves of the Cobra Lily. These data come from Gotelli &amp; Ellison’s text book <a rel="nofollow" href="https://global.oup.com/academic/product/a-primer-of-ecological-statistics-9781605350646?cc=ca&lang=en&" target="_blank"><em>A Primer of Ecologisal Satistics</em></a>. Whether or not a wasp visited a leaf during the census was recorded along with the height of the leaf from the ground. The aim is to test the hypothesis that the probability of leaf visitation increases with leaf height.
</p>
<p>
Let’s jump right in and fit the GLM, a logistic regression model
</p>
<figure>
<pre><code data-lang="r"><span>mod</span><span> </span><span>&lt;-</span><span> </span><span>glm</span><span>(</span><span>lvisited</span><span> </span><span>~</span><span> </span><span>leafHeight</span><span>,</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>wasp</span><span>,</span><span> </span><span>family</span><span> </span><span>=</span><span> </span><span>binomial</span><span>())</span><span>
</span><span>summary</span><span>(</span><span>mod</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Call:
glm(formula = lvisited ~ leafHeight, family = binomial(), data = wasp)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.18274  -0.46820  -0.23897  -0.08519   1.90573  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -7.29295    2.16081  -3.375 0.000738 ***
leafHeight   0.11540    0.03655   3.158 0.001591 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 46.105  on 41  degrees of freedom
Residual deviance: 26.963  on 40  degrees of freedom
AIC: 30.963

Number of Fisher Scoring iterations: 6</code></pre>
</figure>
<p>
Now create a basic plot of the data and estimated model
</p>
<figure>
<pre><code data-lang="r"><span>## some data to predict at: 100 values over the range of leafHeight
</span><span>ndata</span><span> </span><span>&lt;-</span><span> </span><span>with</span><span>(</span><span>wasp</span><span>,</span><span> </span><span>data_frame</span><span>(</span><span>leafHeight</span><span> </span><span>=</span><span> </span><span>seq</span><span>(</span><span>min</span><span>(</span><span>leafHeight</span><span>),</span><span> </span><span>max</span><span>(</span><span>leafHeight</span><span>),</span><span>
</span><span>length</span><span> </span><span>=</span><span> </span><span>100</span><span>)))</span><span>
</span><span>## add the fitted values by predicting from the model for the new data
</span><span>ndata</span><span> </span><span>&lt;-</span><span> </span><span>add_column</span><span>(</span><span>ndata</span><span>,</span><span> </span><span>fit</span><span> </span><span>=</span><span> </span><span>predict</span><span>(</span><span>mod</span><span>,</span><span> </span><span>newdata</span><span> </span><span>=</span><span> </span><span>ndata</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>'response'</span><span>))</span><span>
</span><span>## plot it
</span><span>plt</span><span> </span><span>&lt;-</span><span> </span><span>ggplot</span><span>(</span><span>ndata</span><span>,</span><span> </span><span>aes</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>leafHeight</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>fit</span><span>))</span><span> </span><span>+</span><span>
</span><span>geom_line</span><span>()</span><span> </span><span>+</span><span>
</span><span>geom_rug</span><span>(</span><span>aes</span><span>(</span><span>y</span><span> </span><span>=</span><span> </span><span>visited</span><span>,</span><span> </span><span>colour</span><span> </span><span>=</span><span> </span><span>lvisited</span><span>),</span><span> </span><span>data</span><span> </span><span>=</span><span> </span><span>wasp</span><span>)</span><span> </span><span>+</span><span>
</span><span>scale_colour_discrete</span><span>(</span><span>name</span><span> </span><span>=</span><span> </span><span>'Visited'</span><span>)</span><span> </span><span>+</span><span>
</span><span>labs</span><span>(</span><span>x</span><span> </span><span>=</span><span> </span><span>'Leaf height (cm.)'</span><span>,</span><span> </span><span>y</span><span> </span><span>=</span><span> </span><span>'Probability of visitation'</span><span>)</span><span>
</span><span>plt</span></code></pre>
</figure>
<figure>
<img alt="Estimated probability of visitation as a function of leaf height." src="https://i0.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/confidence-intervals-for-glms-darlingtonia-plot-fit-1.png?w=456&ssl=1" data-recalc-dims="1"><figcaption>
Estimated probability of visitation as a function of leaf height.<br>
</figcaption></figure>
<p>
Next, to illustrate the issue, I’ll create the confidence interval the <em>wrong</em> way
</p>
<figure>
<pre><code data-lang="r"><span>## add standard errors
</span><span>ndata</span><span> </span><span>&lt;-</span><span> </span><span>add_column</span><span>(</span><span>ndata</span><span>,</span><span> </span><span>wrong_se</span><span> </span><span>=</span><span> </span><span>predict</span><span>(</span><span>mod</span><span>,</span><span> </span><span>newdata</span><span> </span><span>=</span><span> </span><span>ndata</span><span>,</span><span> </span><span>type</span><span> </span><span>=</span><span> </span><span>'response'</span><span>,</span><span>
</span><span>se.fit</span><span> </span><span>=</span><span> </span><span>TRUE</span><span>)</span><span>$</span><span>se.fit</span><span>)</span><span>
</span><span>## compute a 95% interval the wrong way
</span><span>ndata</span><span> </span><span>&lt;-</span><span> </span><span>mutate</span><span>(</span><span>ndata</span><span>,</span><span> </span><span>wrong_upr</span><span> </span><span>=</span><span> </span><span>fit</span><span> </span><span>+</span><span> </span><span>(</span><span>2</span><span> </span><span>*</span><span> </span><span>wrong_se</span><span>),</span><span> </span><span>wrong_lwr</span><span> </span><span>=</span><span> </span><span>fit</span><span> </span><span>-</span><span> </span><span>(</span><span>2</span><span> </span><span>*</span><span> </span><span>wrong_se</span><span>))</span></code></pre>
</figure>
<p>
and plot the resulting interval
</p>
<figure>
<pre><code data-lang="r"><span>plt</span><span> </span><span>+</span><span> </span><span>geom_ribbon</span><span>(</span><span>data</span><span> </span><span>=</span><span> </span><span>ndata</span><span>,</span><span> </span><span>aes</span><span>(</span><span>ymin</span><span> </span><span>=</span><span> </span><span>wrong_lwr</span><span>,</span><span> </span><span>ymax</span><span> </span><span>=</span><span> </span><span>wrong_upr</span><span>),</span><span>
</span><span>alpha</span><span> </span><span>=</span><span> </span><span>0.1</span><span>)</span></code></pre>
</figure>
<figure>
<img alt="Estimated probability of visitation as a function of leaf height with an incorrectly-computed 95% confidence interval superimposed. Notice the interval exceeds the probability limits, 0 and 1." src="https://i1.wp.com/www.fromthebottomoftheheap.net/assets/img/posts/confidence-intervals-for-glms-add-wrong-interval-1.png?w=456&ssl=1" data-recalc-dims="1"><figcaption>
Estimated probability of visitation as a function of leaf height with an incorrectly-computed 95% confidence interval superimposed. Notice the interval exceeds the probability limits, 0 and 1.<br>
</figcaption></figure>
<p>
That’s problematic because for significant sections of <code>leafHeight</code> our uncertainty interval breaks the laws of probability.
</p>
<p>
So, when creating confidence intervals we should expect asymmetric confidence intervals that respect the physical limits of the values that the response variable can take. If they don’t, then you’ve probably computed them the wrong way.
</p>
<p>
The previous paragraphs walked through a logical reason why confidence intervals are not symmetric on the response scale. The theory behind adding/subtracting two times the standard error is also derived for models where the response is conditionally Gaussian. It doesn’t really work properly at all when the response is not conditionally distributed Gaussian; you only need to realise that a confidence interval that includes impossible values can’t possibly have the coverage properties claimed because some part of it lies in a space of values that just won’t ever be observed.
</p>
<h3>
Confidence intervals the right way<br>
</h3>
<p>
How do we create correct confidence intervals?
</p>
<p>
A simple solution is to create the interval on the scale of the link function and not the response scale. On the link scale, we’re essentially treating the model as a fancy linear one anyway; we asssume that things are approximately Gaussian here, at least with very large sample sizes. Given that assumption, we can create a confidence interval as the fitted value plus or minuss two times the standard error on the link scale, and the use the inverse of the link function to map the fitted values and the upper and lower limits of the interval back on to the response scale.
</p>
<p>
If you paid attention in your stats classes, you might know that the default link for the Poisson GLM is the logarithm link. You might also know that the inverse of taking logs is exponentiation. You may even know that exponentiation is done in R using the <code>exp()</code> function. But what’s the inverse of the logit function, which was the link used in our model for leaf visitation? Even if you knew what the correct mathematical function was, would you know what R function to use for this? And I defy most readers to know what the inverse of the complementary-log-log link function is, which we could have used instead of the logit link in our model. This problem only gets worse when we start thinking about models that walk and quack like a GLM but aren’t really GLMs in the strict sense, but which use families that are outside the usual suspects of the exponential family of distributions.
</p>
<p>
All is not lost however as there is a little trick that you can use to always get the correct inverse of the link function used in a model. (Well, always is a bit strong; the model needs to follow standard R conventions and accept a <code>family</code> argument and return the <code>family</code> inside the fitted model object.)
</p>
<p>
Typically in R, functions that fit generalized models take a <code>family</code> argument and return a <code>family</code> object that we can extract from the model itself. That <code>family</code> object contains all the information we need to create proper confidence intervals for GLMs and related models.
</p>
<p>
For the logistic regression model we fitted earlier, the family object is the same as that returned by <code>binomial(link = 'logit')</code>, and we can extract it directly from the model using the extractor function <code>family()</code>
</p>
<figure>
<pre><code data-lang="r"><span>fam</span><span> </span><span>&lt;-</span><span> </span><span>family</span><span>(</span><span>mod</span><span>)</span><span>
</span><span>fam</span><span>
</span><span>str</span><span>(</span><span>fam</span><span>)</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">Family: binomial 
Link function: logit 

List of 12
 $ family    : chr &quot;binomial&quot;
 $ link      : chr &quot;logit&quot;
 $ linkfun   :function (mu)  
 $ linkinv   :function (eta)  
 $ variance  :function (mu)  
 $ dev.resids:function (y, mu, wt)  
 $ aic       :function (y, n, mu, wt, dev)  
 $ mu.eta    :function (eta)  
 $ initialize:  expression({  if (NCOL(y) == 1) {  if (is.factor(y))  y &lt;- y != levels(y)[1L]  n &lt;- rep.int(1, nobs)  y[weights =| __truncated__
 $ validmu   :function (mu)  
 $ valideta  :function (eta)  
 $ simulate  :function (object, nsim)  
 - attr(*, &quot;class&quot;)= chr &quot;family&quot;</code></pre>
</figure>
<p>
If you look closely you’ll see a component named <code>linkinv</code> which is indicated to be a function. This is the <em>inverse</em> of the link function. The link function itself is in the <code>linkfun</code> component of the family. If we extract this function and look at it
</p>
<figure>
<pre><code data-lang="r"><span>ilink</span><span> </span><span>&lt;-</span><span> </span><span>fam</span><span>$</span><span>linkinv</span><span>
</span><span>ilink</span></code></pre>
</figure>
<figure>
<pre><code data-lang="text">function (eta) 
.Call(C_logit_linkinv, eta)
</code></pre></figure>