# Generalised Linear Models (GLM)
It is possible to generalise the standard regression models you already know so that you can handle cases where your residuals are not Normally distributed. Remember that *it's the errors that matter*: in previous statistics classes you may have been told that the distribution of the response variable is what matters the most, but this was an over-simplification. These kinds of models are not just more useful because they fit your data well, they also provide predictions that make more sense. Generalized Linear Models (GLMs) can be used to predict integers (*e.g.*, there should be four or five species over there, not 4.5), fractions and percentages (*e.g.*, 99% of people want pizza for dinner, not 101%), and binary variables (*e.g.*, the patient will live or die, there's nothing inbetween). GLMs can be a bit scary at first: focus on the concepts and make sure you've nailed down the two pieces that are at play: the *error structure* and the *link function*. There is a lot of extension material here: I would be quite happy if, for example, you focused only on the Poisson example, or watched the lecture video (which gives a very high-level overview) and then attempted only the exercise and/or to follow along with the code examples below. If you have been following the course so far and found it easy, then great, attempt these exercises, but otherwise consider going back through all the material covered so far, perhaps including the extension exercises there, and reviewing that.
## Link functions
Not all data are Normally distributed. A good example is *count data*—data where you have counted whole numbers of things, like the number of species or individuals in a quadrat. Such data can't be Normally distributed: you can't have 'half a species' or 3.5 individuals, and the continuous values that a Normal distribution produces are inapproriate for such data. Another example is fraction or percentage data. Normal distributions are unbounded, but it's impossible for someone to "give 110%", and so it would be better to model our data as coming from a distribution that was bounded between 0 and 1 (0% and 100%).
*Link functions* give us exactly that power. The equation for a link function looks quite simple (depending on how comfortable you are with Greek letters), but it's important you nail down the concepts:
$$
\eta = g(\mu)
$$
Where $\eta$ is the *linear predictor*, $g$ is the *link function*, and $\mu$ is your prediction (as defined earlier in the course, such as in @eq-regression_equ). The link function is what allows a GLM to make predictions that aren't Normally distributed: it transforms your predictions so that a linear model can give a non-linear response. This transformed prediction ($\eta$; the linear predictor) is then what is fit in your model. Thus if you want to take the parameter coefficients of your model and understand what they mean for your real data, you must take the *inverse of your link function* to back-transform your data.
The simplest link function is the *identity link*—it just leaves your data as-is. Nothing changes with the identity link. The link function we'll be working with today for the Poisson distribution is the *log link*—it's the natural logarithm of your data. It's well-suited to data that are roughly log-normally distributed—like count data. So if you fit a GLM with a *log link*, the parameters would all be estimated in 'log-space', and to make sense of them in your original data you would have to put the data through the inverse of the log-transformation (*i.e.*, take the natural exponent of them, something like $e^y$). If you want to think of the link function as just a fancy way of transforming your response variable, that's fine for now.
## Error structures
Residuals (errors) are not always Normally distributed. In an ideal world, we would fit models that wouldn't predict 'half a species' (0.5) or any kind of non-integer number of species. GLMs let us specify this by defining the *error structure* of our data. It's simply a way of telling `R` what kinds of errors to expect, and it lets `R` adjust its model-fitting procedures accordingly^[For a full understanding of what error structures do, you need to understand the algorithms by which GLMs are numerically fit. To do so is beyond the scope of this course; if you're interested, I suggest reading either 'Extending the Linear Model with `R`' and/or 'Likelihood', the citations for both of which are in the forward. You can also ask me about it, but it will be tricky for us to have a conversation about it if you haven't taken my programming class, where I teach basic numerical methods.].
Error structures go hand-in-hand with link functions, and the choice of error structure is often what determines the link function you use. Putting that a little more formally, while it is possible to mix-and-match link functions and error structures, there is almost always a *canonical link* function associated with each error distribution that you're better-off using. These links have been derived analytically, and are the pairings most likely to give you precise, unbiased predictions. My advice is to *let the kind of data you have determine your error structure*, and then always to use the canonical link function, although if you are interested it is possible to statistically test which link is most appropriate using $F$-tests (ask me if you are interested). The canonical error distribution for the Poisson distribution that were working with today is linear—the variance is directly proportional to the mean. I outline this in more detail below.
## Practical example: count data and the Poisson family
Count data are weird. As we've already discussed, you can't have negative counts (they're *bounded*), and they can only be *integers* (whole numbers, not fractions—*discrete* data). Additionally, the variance of count data tends to increase with its mean: they show *non-constant variance*. Putting that another way, if I asked you to guess how many people were in the biology department right now, you might guess something like $100 \pm 10$. If I asked you to guess how many people in the university right now, you might guess something like $5000 \pm 1000$—you definitely wouldn't guess $5000 \pm 10$ because you'd understand, intuitively, that the error in your estimate has increased with the mean.
The error distribution that has all these properties is called the *Poisson distribution*, and is named after the French mathematician Siméon Denis Poisson. Unlike the Normal distribution, which has two parameters (the mean—$\mu$—and the variance—$\sigma^2$), the Poisson has only one parameter—$\lambda$ ('lambda')—which parameterizes both the mean and variance of the distribution. As a result, it's almost always positively skewed (it has a long positive tail). It also only ever predicts whole numbers (it's discrete), which makes it perfect for count data. You can use something like `hist(rpois(1000, 5))` to plot a histogram of 1000 draws from a Poisson distribution with a $\lambda$ of 5; play around with different $\lambda$ values to get a feel for the distribution.
Let's simulate some data that represents the abundance of a species across an temperature gradient. The species is more abundant when it gets warmer, and its abundance across sites is drawn from a Poisson distribution. First of all, we'll model the data using a standard linear regression, just to see what a poor job it does of explaining this data.
```{r}
# Simulate data
temp.sites <- seq(-25, 25, by=.5)
lambda.sites <- temp.sites + abs(rnorm(length(temp.sites)))
abundance <- numeric(length(lambda.sites))
for(i in seq_along(abundance))
abundance[i] <- rpois(1, lambda.sites[i])
abundance[is.na(abundance)] <- 0
# Plot data
plot(abundance ~ temp.sites)
# Model data (poorly)
bad.model <- lm(abundance ~ temp.sites)
# Plot the model predictions
abline(coef(bad.model), col="red", lwd=3)
# Inspect residuals of our model
hist(residuals(bad.model))
```
Notice how the residuals in our model are very clearly not Normally distributed. OK, now let's fit a better model to our data, compare its model performance, and then plot the predictions from it.
```{r}
# Fit a better model
good.model <- glm(abundance ~ temp.sites, family=poisson)
# Inpsect residuals
hist(residuals(good.model))
# Summarize model
summary(good.model)
```
The residuals look better here than before, but note that our GLM is still getting strange (but not as strange) residuals from all those zero values we have in our model^[If you really want to account for these kinds of effects, you may need to fit a *zero-inflated* model, which I don't cover in this course, or *quasi-likelihood* models, which I do.]. A thing I do want you to notice is that the residuals still look somewhat strangely distributed: by default, `R` plots the same kinds of residuals as it does for a standard regression. If you want to examine the *working* residuals (what the GLM is fit to, once your data have been put through the link function) check out the help file for `residuals.glm`^[There is a good discussion of this in 'Extending the Linear Model in `R`'.]. When you're looking at the model summary, ignore anything related to *deviance* as we're going to cover this in a later session—for now, it's sufficient to think of these as related to the $r^2$ values you know from standard regression. Since our data are not Normally distributed, we can't use $r^2$ values any more, but deviance is a generalization of the same essential concept. Our coefficients are all on the same scale as the link function, so to work with them, we need to transform everything. So let's do that, and compare the fit of our models.
```{r}
# Make model predictions
predictions <- data.frame(temp.sites=seq(-25, 25, .01))
predictions$good <- predict(good.model, predictions, type="response")
predictions$bad <- predict(bad.model, predictions)
# Compare model predictions
plot(abundance ~ temp.sites, pch=20)
with(predictions, lines(good ~ temp.sites, col="red", lwd=3))
with(predictions, lines(bad ~ temp.sites, col="blue", lwd=3))
```
First of all, notice how, when transforming my coefficients, I have to add all of them together and *then* back-transform. You can't take the exponent of the slope, multiply that by the temperature, and then add that to the exponent of the intercept. *The entire regression equation must be transformed together*; this is simply how logarithms work, and if the math is not obvious to you just remember then just follow my rule in italics. Notice how I can get `R` to do my back-transforming for me by specifying `type="response"` when using the `predict` function. Finally, there are two things I hope are clear from the plot of our model predictions: (1) the Poisson GLM is doing a much better job of capturing the shape of our data, and (2) it is curved when plotted in linear space because our link function is a logarithm.
## Reviewing the Binomial distribution
The *Binomial* distribution is one of the most important distributions in mathematics. It underlies a lot of work in the theory of combinations and permutations, which, in turn, underlies a lot of the work in cryptography. You might find it surprising that there are many questions we can ask using special cases of the Binomial distribution that we can't actually find an answer for!
The Binomial distribution gives the probability of a given *event* happening a given number of times out of a given number of *trials*. It can be used to answer questions like "what's the probability I will get five heads if I toss a coin thirty times" or "how many of these seventy patients would we expect to survive given a 20% mortality rate". Thus the Binomial distribution is defined by two parameters: the probability of success, $p$, and the number of trials, $n$. In the special case where only a single trial occurs, the Binomial distribution simplifies down to something called the *Bernoulli* distribution. Note that the Binomial distribution is, itself, a special case of the *Multinomial* distribution; the multinomial is used in cases where there is more than one kind of event that can happen, such as questions like "what's the probability I'll roll a 1, 2, 3, 4, or 6 if I roll a die five times".
The Binomial distribution is well-suited to cases where we have *bounded* data, since it can only produce estimates between 0 and 1 (inclusively). Examples of this kind of data are percentages (99% of people love statistics!) and fractions (49 out of 50 people in the class love statistics!). The Binomial also has non-constant variance: it shows a quadratic relationship with the value of mean (lesser closer to 0 and 1, higher closer to 0.5). This often maps onto how we would want our model to perform: our model probably does worse with intermediate cases, but predicts the extremes much better. @fig-variances shows the expectations for the variances of the three distributions we are covering in this class.
{#fig-variances}
## The logit link function and logistic regression
*Logistic regression* is the term frequently used to describe Generalized Linear Models modeled using the Binomial family. This is because it makes use of the *logistic* equation for modeling a probability of success:
$$
p = \frac{e^{\mu}}{1 + e^{\mu}}
$$ {#eq-logit-nonlinear}
where $p$ is the probability of success, $e$ is Euler's constant, and $\mu$ is our mean response, is non-linear and so cannot be fit using standard linear methods. We can, however, linearize @eq-logit-nonlinear through some fancy algebra^[...which I'm not going to show you, but if you're interested ask...] into the following equation:
$$
ln(\frac{p}{q}) = \mu
$$ {#eq-logit}
Where $p$ is the number of successes, $q$ is the number of failures ($1-p$), and $\mu$ is still the mean response. $ln\frac{p}{q}$ is called a *logit transformation* of $p$, and so @eq-logit is the definition of the *logit link* for the Binomial family. Hopefully you're not a betting person, but if you pretend you are it can make the logit transformation easier to understand: when a bookie says "2 to 1 odds" they mean "$p=\frac{1}{3}$ and $q=\frac{2}{3}$". Confusingly, everyone other than statisticians expresses probability as odds ratios, which are precisely what the logit transformation does. Thus only those who have taken a statistics class are confused by the logit transformation, and everyone who hasn't taken a statistics class grasps it essentially intuitively^[There is perhaps a proof of the reality of quantum statistical physics in this statement.].
You will remember, of course, that there are two parts to a Generalized Linear Model: the link function and the error structure. Thus, to fit a Binomial GLM, we will need to use both a logit link and the error structure outlined in @fig-variances.
## Special cases of the Binomial
Classes in Generalized Linear Models can become very complex when the Binomial family is tackled, because there are two 'special cases' of the Binomial family that are usually treated as separate problems. They're not, and in the hopes of simplifying things for you I want to present them as what they really are: special cases of the application of a Binomial GLM.
**Binary response variables.**
Data are often presented as binary response variables: infected or healthy, a coin came up heads or it didn't, a student liked the class or they didn't. Such binary data can be modeled with a Binomial GLM because such data can still be modeled as drawn from a Binomial distribution. The only difference now is that, in each case, the number of trials ($n$) is 1: that is the only difference. Indeed, in many cases data that could be presented in the format described above (successes vs. failures) can be re-coded as a binary response variable by 'expanding' out the data. For example, if I tossed 100 coins and got 48 heads, that could either be coded as 48 successes and 52 failures (two pieces of data) or expanded out to be 100 pieces of data, each describing the outcome of a single coin toss (success/fail, head/tail, 1/0, or true/false). It would give exactly the same answer in a GLM because it's exactly the same mathematics.
**Percentage response variables.**
Percentages are really just fractions: a number of times something happened divided by the total number of times something could have happened^[Statisticians never give 110%.]. Knowing how many times something could have happened is extremely important. Imagine I told you I tested my die and it rolled a six 50% of the time: would you be surprised? You probably would be if I told you I'd rolled it 1000 times, but if I'd only rolled it twice you probably wouldn't be. This problem permeates through into all statistics: if you don't know the numerator (1 or 500 in my die example) and the denominator (2 or 1000 in my die example) underlying some percentage data, you are missing critical information that will help you do your job better. Thus, whenever possible, *work with the raw counts that generated the percentages*. If you know the percentage and the number of trials, it's possible to reconstruct the raw data ($successes=\frac{n}{\%}$) and you should do so.
We do not, however, live in a perfect world, and the reality is that as a practicing statistician you will rarely be called in to analyze datasets that are perfectly formatted and set up exactly as you would wish. In the event that you are given a percentage, and do not know the underlying counts that it came from, it is still possible to work with that data. However, it's not possible to model those data using a Binomial GLM: without any information about the number of trials ($n$) it simply isn't possible. Your best bet is to transform your response variable, using something like an arc-sine transformation, and model that using a standard linear regression^[I suppose you could put your data through an arc-sine link function within a GLM framework, but that's not really going to give you any additional information.]. If I may give you a word of caution and/or advice, however: if the person who gave your the percentage won't tell you the number of trials, be very hesitant about moving forward analyzing that data. I can't really imagine how that's the kind of data one could lose, and I would want a very good explanation of why that data wasn't available. I am perhaps a tad paranoid because I briefly worked in industry, where such information would be witheld deliberately because of... politics.
## How to tell if a GLM fits well
For multiple regressions, you already know the answer: the $r^2$ of a model is its ability to predict the data on which it was fit. A number closer to 1 is better than a number closer to 0, and what value is absolutely good (is 20% good?) depends on your field and question. Of course, if a model isn't statistically significant (as measured by an $F$-test; see below for a refresher) then it doesn't really matter whether the $r^2$ is good or not, but if you find you have a good $r^2$ and no significance then you've almost certainly done something very peculiar and so should think about whether you have enough data to fit a model full stop.
For everything else, the problem is more complex. There is no such thing as an $r^2$ for anything other than a multiple regression (*i.e.*, a GLM with the identity link). People commonly calculate a form of pseudo-$r^2$ where they take model predictions and correlate the original data with those: that's fine and everything, but it's not an $r^2$ and so it's not clear what that means. Real statisticians use the *deviance* of a model: it is what is minimized to fit your model to data, and its definition is determined by your error distribution. @tbl-deviances lists the most common ones; you don't need to memorize them, but it can be useful to know what they are.
You will have noticed that the deviance for Normal errors is just the sum of the squares; the sums of squares are, indeed, a special case of the deviance. That's why bad statisticians often plop the deviance into the $r^2$ formula I've already given you, and use that to calculate an $r^2$. Such people are silly, because they don't propagate the changes through into the other parts of the $r^2$ equation (because it can't be done) and so real statisticians like to make fun of such people in the pub after conferences. Remember that, during the fitting process, each datapoint is weighted according to the *error distribution* (see @fig-variances): thus each datapoint is no longer equally important, but this change isn't sensibly propagated through the mathematics if you just divide the null by the residual deviance. Thus the ratio of the two numbers is not quite meaningless, but it's not an $r^2$ value. Real statisticians report the null deviance (the amount of deviance in their data to begin with), their residual deviance (the amount left behind by their model), and plot the raw data with the model predictions. This is sufficient.
| Error structure | Deviance equation | Variance function |
| :--- | :--- | :--- |
| Gaussian (Normal) | $\sum(y - \bar{y})^2$ | $\propto 1$ (*i.e., a constant*) |
| Poisson | $2\sum(y log(\frac{y}{\mu}) - (y-\mu))$ | $\propto \mu$ |
| Binomial | $2\sum(y log(\frac{y}{\mu}) - (n-y) log(\frac{n-y}{n-\mu}))$ | $\propto p(1-p)n$ |
: **Some common deviance equations.** The error structure, deviance equation, and relationship between the predicted response and its error, are given in each column, respectively. $y$ is the response variable, $\bar{y}$ the mean value of $y$, $\mu$ the predicted response value, $n$ is number of trials, $p$ is the probability of success, and $log$ is the natural logarithm. Note, in particular, the last column: this expected variation is why it's often not possible to get an $r^2$ from non-Gaussian GLMs: the expectation of error is not the same for all predictions, thus not all error is equally important, thus not all deviance is equally important, thus you cannot simply sum deviance as where the deviance 'comes from' matters. No equation is given for the logistic deviance, since this is a special case of the Binomial. See also @fig-variances. {#tbl-deviances}
For mixed effects models, which we will be learning about at the end of the course, this problem is a lot more complex. We'll cover that when we cover them, but I want to flag now that anyone who quotes you an $r^2$ for a mixed effects models model is either ignorant or being deliberately misleading. Anyone who quotes you two $r^2$ values likely knows what they're doing, but I would still argue you should be cautious in interpreting those values as they don't map onto $r^2$ as it has been defined in this course.
## Overdispersion and quasi-likelihood
Each of the GLM families we've encountered have equations for their variance (see @fig-variances and @tbl-deviances). Sometimes, the left-over variation in your model is going to be much greater than the number of degrees of freedom you've got left-over in your model. In cases such as this, we would describe these data as *overdispersed*: the unexplained variation is much greater than we would expect given our statistical model. Our approach to dealing with this in a standard linear model (*i.e.*, a GLM from the Normal family) is to fit a different error structure—to fit a Poisson or Binomial GLM. So what do we do if a Poisson or Binmoial GLM still isn't up to scratch?
One solution is to fit an even more complicated model, either in terms of its underlying formulation^[Classic ones in ecology are 'zero-inflated' models or mixture models that allow your data to come from more than one statistical distribution. I prefer mixture models, but they're both outside the scope of this course I'm afraid. You can also, in a 'Generalized Least Squares' model, fit more complicated error distributions, or in a mixed effects model soak up the variation with extra random effects. We cover these towards the end of this course.] or by adding additional explanatory variables, but the other is to give our error distributions even more wiggle-room. What we can do, mathematically, is to assume that our data are drawn from more than one Poisson or Binomial distribution, each with their own separate error distribution. In principle, the data can be drawn from an infinite number of these equations (and the mathematics are derived under that assumption), but in practice `R` will try and estimate the number of distributions for us. Because we're sort-of guessing the distributions, and there could be so many of them, we call such error distributions *quasi-likelihood* models. We're assuming that our data are drawn from a sort-of mixture of distributions all centered around one defined distribution, and it is literally because the whole thing is so vague that this is called quasi-likelihood. The approximate number of distributions roughly maps onto something `R` calls the *dispersion parameter*, but I can't go into much more detail about this in this course^[If you're interested, 'Extending the Linear Model in `R`' has quite a nice section on this.].
It's a mathematical trick to make sure our estimates aren't biased by all this weird, unexplained error, and so we flag it as a quasi-likelihood model to let everyone know we know we're doing something a bit odd. Mathematically, this is very similar to the trick we pull when we fit so-called 'mixed' models (see later in the course), and just as with mixed models there's nothing inherently wrong with it but it is still a 'mathematician's trick'. We just have to remember there's a lot of unexplained variation in our model, and so we would ideally like to find some more explanatory variables to soak that up if we could.
There is no set-in-stone rule for when you need to fit a quasi-likelihood model, but you should (in my opinion) definitely fit one if you have a greater residual deviance than your residual degrees of freedom. I cannot emphasize enough, however, that the best solution is to find more explanatory variables that will soak up the variation in your response variable. There must be *some* reason that your data are the way they are, and while it's always possible to come up with a statistical fix it's much more satisfying to figure outhe the *biology* underlying your problem.
## Hands-on with Binomial regression
OK, let's go fit a model! We're going to simulate some data about whether students are happy or sad depending on the time-of-day of the class they are sat in. First of all, here's some code to simulate that data.
```{r}
# Inverse logit function for simulation
inv.logit <- function(x) exp(x)/(exp(x) + 1)
# Setup some times at which to measure mood
time <- seq(-3,3,.2)
# Draw successes (happy) from distribution
happy <- rbinom(length(time), 30, inv.logit(time))
# Sad students aren't happy
sad <- 30 - happy
# Merge that data into a two-column explanatory variables
mood <- cbind(happy, sad)
```
Once again, you really don't need to worry about the details of the simulation here, but you do need to take note of the last line. A Binomial GLM requires information on how many successes (happy students) there have been, as well as the number of failures. Thus we use the function `cbind` to bind together two variables and make each a column in a new variable. Once we have this, it's remarkably simple to fit a Binomial GLM:
```{r}
model <- glm(mood ~ time, family=binomial)
summary(model)
```
The code is otherwise exactly the same as the Poisson example we saw last time. For all the reasons described above, we can't really compute an $r^2$ for this model. We can, however, get a fairly good impression of how well it's doing by examining the `Null deviance` and the `Residual deviance`. Your simulated data will likely be slightly different, but in my dataset the null deviance is roughly 463 and the residual deviance is roughly 22. Thus we've explaining something like 95% of the deviance in this data ($\frac{null-residual}{null} \approx \frac{SST-SSE}{SST} = \frac{SSM}{SST}$).
For completeness, and because this is a fairly simple model, let's get a significance test on the differences between the variances of this Binomial model and a null model with no terms. In `R` the symbol for "only fit a model with an overall mean" is `1`, which you can remember because the term `(Intercept)` is always the first term in the model output.
```{r}
null <- glm(mood ~ 1, family=binomial)
anova(model, null, test="F")
```
This is the $F$-test that was described in @eq-aov_equ, and is the same code as in @sec-anova_practical. The only difference is we have to explicitly tell `R` to perform an $F$-test here or it won't calculate a $p$-value for us. Were we working with a Poisson GLM, we would have to tell `R` to perform a `Chi` ($\chi^2$) test here. The $F$ distribution is essentially a normalization of the $\chi^2$-distribution, and as the number of degrees of freedom in $SSM$ approaches infinity the two converge. Thus there's not really a tremendously good reason that you *need* to switch between the two when you're testing significance in GLM, but `R` will shout at you if you don't and it's always good to avoid being shouted at!
I want to remind you that I only really think nested ANOVAs of this type should be used when you have very specific hypotheses you want to test. In standard statistics classes, we would be constantly performing these tests to get the significance ($p$-value; see session 1) of each parameter/term in our model. I'm very specifically asking you not to do that here: please just trust me, and wait until the next session when we will go through, in mind-numbing detail, how to figure out when terms in your model are significant or not.
Finally, there is very little evidence that our model is overdispersed (the residual deviance is vastly smaller than its degrees of freedom), but if there were then we would fit an overdispersed model like this:
```{r}
quasi.model <- glm(mood ~ time, family=quasibinomial)
summary(quasi.model)
```
The interpretation of the model would be exactly the same as before. The only differences, if you squint at the output very carefully, is we no longer have a definition of the model's AIC and we have an estimate of the dispersion parameter I mentioned above (greater values mean more of a quasi-correction applied, and 1 is the baseline in a standard likelihood model). We don't know what AIC is yet, but when we do, this will become important to know.
## Exercises
```{r}
# Load the GLM Poisson data in like this
data <- read.csv("ants.csv", as.is=TRUE)
# Plot the data out like this
with(data, plot(species ~ elevation))
# Wouldn't it be fun to log the vertical axis?
with(data, plot(species ~ elevation, log="y"))
# Load the GLM Binomial data in like this
data <- read.csv("crabs_simplified_logistic.csv", as.is=TRUE)
# Record which crabs are damaged like this
data$damaged <- with(data, as.numeric(ifelse(legs != 8 | pincers != 2, TRUE, FALSE)))
# Plot the data out like this
with(data, boxplot(diameter ~ damaged, horizontal=TRUE))
# ...or like this - note the use of jitter on the response variable
with(data, plot(jitter(damaged) ~ diameter))
```
1. Before doing anything else, step through all of the code above and make sure you understand what it is doing. Run all the practical examples and make sure you can follow through, from start-to-finish, what every step of the analysis is doing.
2. (Extension) Our GLM Binomial data comes from one of my undergraduate projects looking at the impact of tourist crab-catching on crabs^[Pearse *et al.* 2014; arXiv:1404.0290.]. There is no need to thank me for giving you this dataset to analyze twice: I know how excited you must be to see it again, and the look of joy on your faces is thanks enough for me. Each row is a crab that was caught and had its carapace `diameter` measured, along with its sex and how many limbs it had. The code above creates a new variable^[Gosh, isn't `R` exciting!] that reports whether a crab is missing a limb and so is damaged^[If you are interested, ask me about the time I had to give the police a report on the welfare of crabs as a result of this study. It's a tale filled with drama—mostly because it involved a play.].
a. Fit a regular regression to these data, using `damaged` as the response variable and `sex` and `diameter` as the explanatory variables. Inspect the histogram of the residuals. Does this model fit the data well? How can you tell?
b. Fit a Generalized Linear Model with a Binomial family to this data. Describe, in a few sentences, what the data show. Hint: these data are not grouped into success/failure like the example data, so it's much easier to fit this model—just use `damaged` as a response variable and don't mess around with `cbind` and the like.
c. Plot the data and a smoothed prediction from your model. Hint: if you used `boxplot` to make your plot, you will want to add 1 to the intercept of your smoothed line as `R` plots boxplots starting at 1 on the y-axis. This sounds worse than it is: you literally do something like `lines(pred$diameter, pred$damaged+1)`^[Incidentally, as far as I know this kind of plot is impossible to create in `ggplot2`. Did I mention that `ggplot2` was written by Hadley Wickham and that I am being 100% completely serious when I say that he is absolutely my favourite programmer in the world and I think he really thinks through everything he does very carefully?]. I mention this solely because many people quite like the look of a boxplot for these kinds of data.
3. (Extension) Our Poisson GLM data come from Gotelli & Ellison (2002; Ecology 83 1604–1609). They describe the number of ant species found in a series of surveys in different habitats. Your task is to analyse what the drivers of diversity are in these data.
a. Fit a regular regression (*i.e.*, with an identity link) to these data with only additive terms for the explanatory variables `latitude`, `elevation`, and `habitat`.
b. Fit a Generalized Linear Model with a Poisson family to this data. Describe, in a few sentences, what the data show.
c. Plot the data and a smoothed prediction from your model. Hint: above I give you code, using the `predict` function, that will give you predictions from a model across a single explanatory variable. To get this to work in a model with two explanatory variables, you will need to give two explanatory variables, only one of which will vary. So something like `data.frame(exp1=seq(-10,10,by=.1), exp2=10)` will give you what you want.
