# Multiple models
This is the third session of this course and, in it, I'm going to show you how to extend the regression/ANOVA models you already know to handle more than one explanatory variable. Conceptually, this isn't much of an extension beyond what you already know, and so it should form a sort-of review of the material from the previous section. It will also provide us with an opportunity to consider model and variable, in which we choose from among a (potentially dizzying) array of possible explanatory variables and statistical model forms. I will end with an overview of extension material—generalised linear models—that I would encourage you to at least glance through because it is advanced material that is normally too complex to quickly cover but, because of the way this course is structured, may seem remarkably straightforward to you now.
## Multiple regression
### ANCOVA? I 'ardly knew 'er!
Multiple regression is often called an ANalysis of COVAriance (ANCOVA). I personally think it's a tremendously unhelpful term^[For reasons I won't go into, because then I'd have to explain what it is, which would confuse you...], but I'm mentioning it here in case you encounter it in some other textbook. A *multiple regression* relies on exactly the same underlying mathematics as a standard regression/ANOVA, but now there can be more than one explanatory variable. Critically, it is not (as some call it) a '*multivariate regression*'; if multivariate regression is anything, it is regression with two response variables, not two explanatory variables, and there is no grounds for calling a multiple regression a multivariate regression that doesn't also apply to linear regression^[It feels good to get all that off my chest.].
Before I give you a general definition of a multiple regression, let's look at an example with one discrete/categorical variable and another continuous variable.
$$
\mu = a_{group} + bx
$$ {#eq-ancova}
Where the terms are exactly the same as they were in @eq-regression_equ and @eq-anova_equ: $\mu$ is the model's *predictions* of the response variable, $b$ is the slope of the equation, $x$ is the explanatory variable, and $a_{group}$ is a different intercept for each group of our response variable^[Those of you who are paying attention may have noticed that this model estimates values for each group, not an intercept and groups' differences from that intercept. Thus isn't exactly the model that is fit by default in `R`. This distinction is important, but for the time being please ignore it until we spend a lecture on contrasts.]. Such a model could give us the kind of model we see in @fig-multiple_additive a, and we could calculate its $SST$ in exactly the same way as we did for a regular regression/ANOVA as shown in @fig-multiple_additive b.
::: {#fig-multiple_additive layout-ncol=2}
{#fig-multiple_additive_a}
{#fig-multiple_additive_b}
**A graphical overview of a multiple regression with group-level intercepts.** In (a), the data are shown as black dots, the error ($x - \mu$) in red, and the model predictions ($\mu$) in blue. Each of the three groups of data ('a', 'b', and 'c') has a separate intercept term (described in the text as $y_{group}$), but they all share a common slope. In (b), we see the same data, but now their differences from the overall mean ($\bar{y}$) are represented in orange. There is nothing mysterious about these figures: compare them with figures from the last session to see how they are natural extensions of the simpler regression/ANOVA with which you are already familiar.
:::
What's notable about @eq-ancova is it only allows our *intercepts* to vary, not our slopes (indeed, there is only one slope). But it's possible to imagine a model where only the slopes varied among the groups, and what its $SSE$ would look like, as you can see in @fig-multiple_interactions a. Indeed, you can even imagine a multiple regression model where we allowed an *interaction* between our two explanatory variables: where we allowed each group to have its own slope and intercept. Such an *interaction model* would look something like @fig-multiple_interactions b. The equation for such an interaction model would be quite straightforward, given what we already know:
$$
\mu = a + b_{group}x + a_{group}
$$
Where the only new thing is the term $b_{group}$, which is a different slope parameter for $x$ for each group in our model. Using this framework we can calculate $SSM$ using @eq-SST just as we could for our linear regression/ANOVA models. We'd take our estimates for $SSE$ and $SST$ (see @fig-multiple_additive) in exactly the same way as before. There's no way, without teaching you some linear algebra, that I can prove this to you, but you can probably guess that the situation is much the same for even more complicated multiple regression models that have multiple discrete and multiple continuous explanatory variables.
::: {#fig-multiple_interactions layout-ncol=2}
{#fig-multiple_slopes}
{#fig-multiple_interaction}
**A graphical overview of more multiple regressions with differing slopes and a common intercept, and full interactions (differing slopes and intercepts).** In all sub-figures, data are shown in black, model predictions in blue, and errors (residuals) in red. In (a), we see a regression that is directly comparable with that in @fig-multiple_additive, only now, instead of a common slope and differing intercepts, there is a common intercept and varying slopes. In (b), we see a full 'interaction' model: there are differing slopes *and* intercepts for each group. As described in the text, however, the calculation of $SSE$, $SSM$, and all the rest is exactly the same.
:::
### Which model is best?
Once you've been given the freedom to fit any kind of model you want, you naturally want to know which model fits your data best. Imagine you were trying to model the number of species in study sites as a function of the temperature, elevation, and habitat type (forest vs. grassland) at those sites. We would prefer the simplest model that explains the most variance in our response variable: is that a model with all three explanatory variables, or just temperature? Do we need an interaction term between habitat and temperature, or can we fit a purely *additive* (*i.e.*, no interaction terms) model?
As I hinted in our last lecture, we can use $F$-statistics to test different models and see whether there is a statistically significant difference in the variance explained between models with different explanatory variables. The protocol for this is straightforward but fiddly: compare the extra variance explained by your more complicated model with the left-over error, account for the number of extra parameters you fit, and then test to see whether the difference is significant. In mathematical symbols:
$$
F = \frac{(SSM_a - SSM_n) \times (n-k_a-1)}{SSE \times (k_a - k_n)}
$$ {#eq-aov_equ}
Where $SSM_a$ and $SSM_n$ are the variances explained (the $SSM$s) by the more complicated and simpler models, respectively, and $k_a$ and $k_n$ are the number of parameters in the more complicated and simpler models, respectively. You can then test the significance of the resulting $F$-statistic in just the same way as you did with @eq-f_stat, using $(k_a-k_n)$ and $(n-k_a-1)$ as your degrees of freedom. Indeed, this equation and @eq-f_stat should look similar because *they're exactly the same*: the only reason this 'new' equation looks more complicated is we're now comparing it with a model with more than one parameter. Indeed, we can use this approach to test any *alternative hypothesis* versus any *null hypothesis*, and that's why I've labeled the sums of squares $SSM_a$ and $SSM_n$.
There are two potential pitfalls with this kind of approach: *order matters* when fitting models and variables, and it's easy to suffer from *multiple testing problems*. We'll cover order in our practical example, but we must discuss testing problems now. Our $\alpha_{crit}$ value that we pick makes our analysis sensitive to the number of tests we perform: if we perform 20 $F$ tests at $\alpha_{0.05}$ then we would expect, by chance, that one of our tests would be significant ($\frac{20}{5\%} = \frac{20}{0.05} = 1$). The old-school solution to this problem was to set a more stringent $\alpha_{crit}$ when testing interactions (see below in the practical example), and to a certain extent there's nothing wrong with this. It's not quite that simple to fix, however, because if we set our $\alpha_{crit}$ to be more stringent (*e.g.*, $\alpha_{0.001}$) then we decrease our chances of finding a real relationship (see lecture 1). In the next section, however, I'm going to show you a much easier way to deal with this problem.
### Multiple-testing corrections and *post-hoc* tests
From the argument above, we can see that multiple testing is areal problem. It is essentially obviated by the problems outlined in the next section, which allow us to focus on effect sizes and so move beyond significance testing entirely, but there are situations in which such approaches are not possible and you need to know how to address them. Luckily, the simplest solution to this problem—the *Bonferroni correction*^[Which was not devised by Bonferroni but rather by Olive Jean Dunn who derived everything by applying smoe maths developed by Bonferroni to solve this problem. I cannot help but find myself wondering if the method would be differently named if Professor Dunn hadn't also happened to be a woman and Bonferroni a man, but I'll leave thinking about that as an exercise to the reader...]—is simply to divide the $\alpha_{crit}$ by the number of tests you are conducting. Thus if you are conducting 100 tests and $\alpha_{crit}=5\%$, then your corrected $\alpha_{crit}=0.05\%$.
This is a very simple approach, but also has the potential to throw the baby out with the bathwater in that you end up with very little statistical *power* (called $\beta$ in statistics to match $\alpha$). In other words, $\alpha_{crit}$ becomes so small that you rarely find anything that has a surprising effect. So people have subsequently corrected the Bonferroni correction^[Perhaps 'improved' would be a better term...], either by being very specific about the order in which you conduct tests and using that to help with the problem (the Holm-Bonferroni correction) or using slightly different multipliers (the Šidák correction). These do, indeed, somewhat improve statistical power, but none of them can get around the fundamental trade-off between making a mistake of omission (power) versus comission (spurious detection), and thus are more susceptible to the kinds of multiple-testing problem we are trying to avoid. I would argue, personally, that these corrections are still better than a naïve application of the Bonferroni correction, so I am still in favour of using them if I have to. But, ultimately, it's best to focus on effect sizes and conducting fewer tests to begin with, in my opinion, which we will cover later in this section and so I won't dwell any more on multiple-testing.
### A practical example {#sec-anova_practical}
Following the example I started earlier, we're going to simulate some data that represents an ecological survey. Imagine you've wandered out into two habitats (grassland and forest) and recorded the number of species found in plots of different temperatures and humidities. Don't worry if the simulation code seems complex to you; I'm happy for you to treat it as a black box and focus on the statistics.
```{r}
# Setup our sites
habitat <- rep(c("grassland","forest"), each=50)
temp <- rep(seq(20,40,length.out=50),2)
humid <- rep(seq(0,20,length.out=50),2)
humid <- humid + rnorm(100,sd=4)
# ...adding extra rnorm or humidity correlates too much with temperature
# Simulate some diversity
species <- ifelse(habitat=="grassland",
temp * 2 - humid, humid*4 - temp
)
#...add some variation, and ensure we don't have 'half' species
species <- round(species + rnorm(100, sd=2))
# Plot it out
plot(species ~ temp, pch=20, col=ifelse(habitat=="grassland","red","black"))
plot(species ~ humid, pch=20, col=ifelse(habitat=="grassland","red","black"))
```
If we were to start with a simpler model and move up from there, which would be the *wrong thing to do*, we might fit a model like this:
```{r}
additive.model <- lm(species ~ temp + humid + habitat)
summary(additive.model)
```
If you have read the extension section from the last session, you know what all the $r^2$ terms above mean, but otherwise please don't give them any thought. What might be surprising is that temperature, which we know is a major driver in this data (because we simulated it that way) is not detected as a significant variable in this model. The reason for this is simple: we've fit too simple a model, and it doesn't account for the *interaction* that we simulated in the data. So let's fit a model with some interaction terms in there.
```{r}
complex.model <- lm(species ~ temp * humid * habitat)
summary.aov(complex.model)
```
This model contains a three-way interaction term: the interactions across temperature, humidity, and habitat, and it's represented in the final line `humid:temp:habitat`. If we were writing out the equation for this, it would represent a different slope of $humid \times temp$ for each habitat type. Let's unpick this for a moment: an interaction between humidity and temperature (let's call it $b_{th}$) is multiplied by the product of humidity and temperature: something like $b_{ht} \times humid \times temp$. In this model, there's one of these for each habitat type in the data. Because we used the function `summary.aov`, we're not getting what would otherwise be extremely confusing model output with lots of coefficients: instead we're getting the *deletion p-values* for each term in the model—what we would get from running each term through @eq-aov_equ. You could carry that out, and fit the simpler model, using something like the following:
```{r}
# Notice the neat "^2" trick
simpler.model <- lm(species ~ (temp + humid + habitat)^2)
# ...it's the same as typing...
simpler.model <- lm(species ~ temp*humid + temp*habitat + humid*habitat)
# Test the models (now you see what "anova" is for)
# - does the p-value look familiar to you?...
anova(complex.model, simpler.model)
# We can carry on the exercise...
summary.aov(simpler.model)
```
Back in the 'good ol' days', this is the part where I had to explain to you a process called *stepwise model simplification*. You would slowly go through your model, starting with the most complicated terms (*e.g.*, third-order interactions, like the one we just discussed, then second-order terms that are 'two-way' interactions), removing the 'least significant' (greatest $p$-value) term, making a new model, then removing its 'least significant' term, and continuing until everything was significant. Critically, you would never remove a 'non-significant' lower-order term if a higher-order term depended on it: for example, a significant interaction between `humid` and `temp` would mean that you couldn't remove the additive `humid` or `temp` terms. This is because a more complicated interaction relies upon the additive terms: without them there is no basis for the interaction^[It is difficult to 'prove' this to you without showing you the math. You can get a feel for the philosophy by trying to think through how you would meaningfully interpret an interaction between two things when neither of them exists in a model]. The whole process was extremely boring, and resulted in something called the *minimum adequate model*.
Luckily, we don't live in the 'good ol' days'^[It was never clear to me that those days were particularly pleasant], and so I'm not going to force you to do this. It's important that you understand the principle, however, because in the case that you have only a handful of predictors it's a useful way to explore data. It's also important to be able to set up explicit hypotheses ("does humidity matter?") and be able to test them, and this kind of framework will let you do that. Of course, you must be careful to start with *as complicated a model as you can first*, and work down from that. As you've probably noticed, this can be very tedious to do in practice, and so later in the course I'm going to teach you how to use *scaled variables* and *model averaging* to quickly figure out what's going on in your data.
Finally, a quick caveat: all the statistical models I'm covering in this series are correlative, and so all the math underlying our models assumes that none of our explanatory variables are correlated with one-another. If they are correlated with one-another (and they almost certainly are a little bit) then it can be difficult to unpick which explanatory variable the response variable is more correlated with. This is why, when you type models into `R`, you often find that the *order* in which you type them into your model matters. Thus the following models don't necessarily give the same output:
```{r}
one <- lm(species ~ temp * humid * habitat)
two <- lm(species ~ habitat * humid * temp)
```
This is a consequence of how the numerical approximations in `R` have been programmed to work: they fit the first terms in the model first, and so they will often get 'more' of the variance associated with them than the others. This is only really a concern in quite pathological models, and so you really shouldn't worry about it other than to make sure that the explanatory variables going into your model aren't very strongly correlated with one-another. If they are, then I suggest you look into two options: *Principal Components Analysis* (PCA) to make fewer variables (that aren't correlated) from your variables, and *lasso regression*. But let's save discussion of those methods for another time.
## Variable and/or model selection
We now must tackle the most obvious, and yet somehow the most difficult to answer, question in statistics: what are the variables that matter the most, and what is my best model? We're going to focus on learning three concepts: *stepwise model simplification* (the wrong thing to do), *information theoretic model selection*, and *standard effect size*. It's tremendously important that you pay close attention in this class because the concepts we're covering are quite new in biology^[They have been standard in statistics for decades, so why this is the case is beyond me.] and so you will frequently encounter people who encourage you to do the wrong thing when selecting models. This likely includes your advisors, so I'm going to try and give you the skills to have an informed conversation with them about what is best for you to do with your data. Do bear in mind that, just because I'm only presenting two 'good' ways to make decisions as to what's going on in your data, that doesn't mean these are the only (or, indeed, the best) ways of doing things^[My personal preference, now that you ask, is to fit a single model that contains all possible sub-models, and then use relative variable importance to determine what's going on. This often means fitting a fully hierarchical model in a Bayesian framework, such that the posterior distribution tells me what matters in my data. We'll cover this a little in the next lecture.]^[One other, common, method is to calculate *relative variable importance*: I show you one way, below, to do that, another is an extension of the variance partitioning methods you'll learn in the mixed effects section of this class, and another I won't cover is described in the `R` package `relaimpo`.]. But before we get there, of course, I'm going to remind you of the method you already know—which just so happens to also be the best method there is.
### Sensibly chosen hypothesis tests—the best way
Throughout this class, I've been emphasizing to you that everything in statistics is really just a test statistic. That is, in fact, the case with $F$-statistics. Go back and look at any equation for an $F$-statistic (@eq-aov_equ is probably the best). Notice anything? That's right: it's just the observation ($SSM_a$), minus the expectation ($SSM_n$), divided by the uncertainty ($SSE$), and (just as we have to compare a $t$-statistic to a distribution) then it's multiplied by a correction factor for degrees of freedom. *The best way to test a hypothesis is to fit a simple test statistic comparing your observation to some expectation*, and if you properly define your observation ($SSM_a$) to a decent expectation ($SSM_n$) using an $F$-test you will never go wrong. That's not controversial, and your advisor will be ecstatic if you can present them with a single statistic to test a single hypothesis that you have.
It really doesn't have to be any more complicated than that. The reason I know so many statistical tests is because, ultimately, I'm not a tremendously good biologist! I can so rarely find a way to simplify my question down to a single test statistic that I compensate by fitting ever-more-complicated statistical models. Trust me: keep it simple. Use the `anova` function to compare two models, and have thought enough about the biology that whatever answer you find is interesting. It'll make your papers easier to read, faster to write, and make you feel tremendously smug in lab meetings.
### Stepwise model simplification—the wrong thing to do
Above, and throughout this course, we've discussed how it's possible to compare different models using an $F$-test. This can be used to determine whether there is a statistically significant difference between a more complex model and a simpler one. There is, in principle, nothing wrong with using a test in this manner, and that logic underlies something you may encounter called the *Likelihood Ratio Test*^[The LRT is is 'simply' the ratio of the likelihood of one model to another, and the $F$-test is a special case of the LRT.]. This process is the foundation for what was commonly known in statistics as *stepwise model simplification*.
In stepwise model simplification, you start with the most complicated model possible and then simplify it as much as possible. You find the 'least significant' term in your model (a $p$ value of $0.3$ would be 'less significant' than a $p$-value of $0.2$), drop that term, and then see if there's a statistically significant difference between the simpler model and the original model. You continue this process until you can drop no more terms, and you are then left with something called the *Minimum Adequate Model*. There are ways of doing this in `R`, but none of them are automated, and it is both extremely laborious to do *and extremely wrong*. Here is an (in-exhaustive) list of reasons it is wrong:
* The 'relative significance' of terms in a model is affected by the order in which they were included in the model in ways that bias the process (*i.e.*, `lm(y ~ x + z)` can have slightly different coefficients from `lm(y ~ z + x)`). This is because regressions are based on correlations, thus we cannot infer causation and so two even slightly correlated variances will affect one-another in ways we cannot control for.
* It leads to biased parameter estimates through 'model selection bias': the term is either in the model or it isn't, and so the process conflates parameter (un)certainty with parameter values. Such model uncertainty is addressed in information-theoretic-based approaches.
* The process is rampantly open to multiple-testing bias. This is such a problem that most textbooks tell the reader to use different significance thresholds for different phases in the simplification process. As well as being totally arbitrary (at least the corrections I mention above are grounded in sound statistical theory, warts and all), this breaks down the mathematical underpinnings of the process, and so loses the guarantee that you will always end up with the 'best' model.
* Most statisticians reject the idea that it is a reasonable question to ask whether there is a single 'best' model (see below). Thus, even if there is a 'best' model, because the question was meaningless the answer is too.
I could go on^[The problem is so pervasive, and so blindingly obvious to statisticians (but not the rest of us mere mortals), that a paper entitled "*Why do we still use stepwise modelling in ecology and behaviour?*" (Wittingham *et al.* 2006; J. Anim. Ecol.) was published and has accrued nearly a thousand citations.]. I occasionally encounter people who use the function `step`, which relies on AIC (see below) to conduct stepwise model simplification. This is actually *worse* than the $F$-test based method I outline above, as it has *literally no theoretical underpinnings*. It's an easy mistake to make, but just because it involves calculating AIC doesn't make it an information-theoretic approach.
### Model selection using information theoretic criteria—a good thing to do
Not all models are equally good, but some are pretty much as good as each other. If two models are essentially indistinguishable ($r^2$ of $0.6$ and $0.599999$, for example, if you've read through my extension material), perhaps that reflects *model uncertainty*: you're pretty confident these two models are doing well, but you're not sure which is the best. It's this idea that underpins *multi-model inference*, which is just a fancy way of saying that's we'll use the output from multiple models, each weighted by our relative certainty in each. I'm going to brush over a few details of how this works, but if you're interested there is a wonderfully concise article by Burnham and Anderson ("*Multimodel inference: understanding AIC and BIC in model selection*"; 2004; Sociological methods & research) that is much quicker to read than their 2002 book (that is definitive) on the topic.
The basic principle is we fit all *biologically plausible* candidate models to data, calculate their relative fit to the data, subset our models down to those who do particularly well, and then our results are the averages across that subset of models (weighted according to relative model performance). In much the same way that we might want to understand the error associated with a parameter estimate (see the next lecture for more info on that), this approach lets us deal with uncertainty in model structure (the explanatory variables we put in our model). Our measure of model performance is called an *information criterion*, and the most popular is the *Akaike Information Criterion* (AIC). It's defined as follows:
$$
AIC = 2 \times k - 2 \times log(likelihood)
$$
Where $k$ is the number of parameters in your model, $log$ is the natural logarithm, and $likelihood$ is your model's likelihood (the thing we maximize in statistics; see lecture 1). Like likelihood values, AIC values are all relative: the absolute AIC value has no meaning, and you can't compare AIC values across different response variables. The number $2$ means that, essentially, models with the same number of parameters that are within about 2 AIC units of each other are often considered by people to be roughly equivalent^[The origins of the $2$ in the equation are... historical, but it's notable that most statisticians like things that are within two orders of magnitude of the maximum likelihood estimate (MLE). This is, in fact, why we tend to use an $\alpha_{crit}$ of 0.05; that corresponds to two log-units from the MLE.]. We call such differences in AIC values $\delta AIC$: the $\delta$ stands for 'difference in'. There is a correction factor for AIC if you are worried you have too little data (which produces $AIC_C$), and a similar metric derived under (but not restricted to use with) Bayesian principles ($BIC$). This is not an esoteric metric; the original paper describing it has been cited over $38,000$ times, and Akaike won the mathematician's equivalent of the Nobel Prize for it^[Well, I'm not certain, but he won so many prizes I find it hard to believe this isn't true.].
It's a lot easier than it used to be to make use of AIC-based model selection. Below I give you a brief example of how to do it. My main pieces of advice would be:
1. `dredge` fits all the possible subsets of your 'top model', and is named to be derogatory, so make sure all the models you're including make sense. The whole process is based on the idea that you check *plausible* models.
2. Trust the *effect size* and *importance* (% of models a variable was included in, weighted by the relative AIC of the models) more than the pseudo-$p$-values you get out
3. Use the *full* estimates, not the estimates calculated only for models where a variable is found. This will *shrink* the value of the estimates such that variables that are only in a few models will not have large effect sizes
4. Scale your variables (see below) before using this. Then you can compare the effect sizes in terms of absolute importance.
5. Check what models come out best. Sometimes all the models are variants on a theme, and sometimes there are multiple themes. Averaging across two totally different kinds of models will give you weird answers: if there are two equally likely kinds of explanation for what's going on, this is exciting and means your paper will have a greater contribution to the field! Don't brush away the detail!
```{r}
# Load package (named for Moomin cartoons)
library(MuMIn)
# Simulate some data (note for later: they have a mean of 0 and SD of 1)
a <- rnorm(100); b <- rnorm(100); c <- rnorm(100); d <- rnorm(100)
y <- a - 0.5*b + 2.7 + rnorm(100)
# Fit our top model, then 'dredge' all the subsets
top.model <- lm(y ~ (a + b + c + d)^2, na.action=na.pass) # what does ^2 do?...
models <- dredge(top.model)
# Subset our models to contain only those within 4 AIC units of the best model
# - this is the standard number, but there's no reason to use one over another really
# - then summarize across the average of these models
summary(model.avg(models, subset=delta<4, fit=TRUE))
# Get the AIC importance values (the 'summed weights')
sw(model.avg(models, subset=delta<4, fit=TRUE))
```
How interesting: `a` and `b` come out on top, but we get some (spurious) associations with other variables^[Note that this problem of spurious interactions driving apparent pattern is a trap that stepwise model simplification would fall into. If you wish, see what would happen if you used stepwise model simplification on this example...]. Which spurious associations shine out will vary a little for you (this is a random simulation, after all), but if it were `c` I would read this as "`c` might be important, but much less so than `a` and `b` on the basis of its low estimates in the averaged model". The large number of interaction terms are what throws our importance estimates off^[Cade (2015) Ecology, 96(9): 2370–2382 discusses this. This is a frustrating article to read because the problems he outlines are real, and ecologists were making these mistakes *despite* their solution already being well-known. The solution is to standardize your variables (see next section).], although our coefficient estimates remain true to the relative importance of everything. That the coefficient estimates are appropriate leads us nicely into standard effect sizes...
### Standard effect sizes—a good thing to do
The general public likes to compare the relative importance of variables in terms of the *size* of their *effect*. Things like "smoking a cigar is twenty times as likely to cause cancer than smoking a single cigarette^[I've no idea if this is accurate but I bet it woke you up.]" are intuitive to understand, but, remarkably, it can be tricky to make statements like that. The reason is variables have different *scales* (centers, means) and *variances* (standard deviations, variability). These differences in scale and variance make parameter coefficients impossible to directly compare, and so you can't make quantitative statements about how much more important one is than the other. The following example should make this clear.
```{r}
x <- rnorm(100, mean=0, sd=1)
z <- rnorm(100, mean=10, sd=20)
y <- 3*x + z + rnorm(100, sd=.25)
model <- lm(y ~ x + z)
coef(model)
# Let's do a bit of R-fu to make the plots appear alongside each other
par(mfrow=c(1,2))
plot(y ~ x)
plot(y ~ z)
```
Looking at these plots, you would probably agree that $z$ shows the stronger relationship with $y$, and so I would say that $z$ is more important than $x$ in this model. Yes, the coefficient for $x$ is the correct one (the equation says that $y = 3 \times x + 1 \times z$), but because $z$ has both more variation and a larger overall value, its effect dominates the value of $y$. This intuitive property of importance is, rightly, missing from this model. So how do we put it back in? Simple: we *z-transform* our data to *standardize* it so that each explanatory variables has the same mean ($0$) and variance ($1$; and so a standard deviation of $1$ too). The equation for this standardization is quite simple, and it's important enough that you should learn it:
$$
z(x) = \frac{x - \bar{x}}{SD_x}
$$ {#eq-standardize}
In other words: subtract the mean ($\bar{x}$) and divide by the standard deviation ($SD_x$). Doing this is `R` looks a little weird, but is quite straightforward:
```{r}
s.x <- as.numeric(scale(x))
s.z <- as.numeric(scale(z))
s.model <- lm(y ~ s.x + s.z)
coef(s.model)
```
It is now possible to make a statement like "$z$ is roughly $7$ times more important than $x$" ($\frac{22}{3} = 7$) with confidence. This is because you've standardized each variable according to its natural degree of variation, and so you can say that, within the context of that variation, this is how much each variable matters when it changes. A common example where this matters is time series analysis: students are often confused that they get a very small coefficient for 'year', because their year is a very large number in comparison with their other numbers. When they transform them (*e.g.*, $slope \times 2019$ becomes $slope \times 1$), the values becomes directly comparable and so things begin to make a little more sense. Some people would prefer to use the phrase "*$z$ has a roughly $7$ times greater effect on $y$ than $x$*"; I think what matters more is being clear, and I have happily used both phrasings depending on context. There is a nice discussion on this in Andrew Gelman's "Bayesian Data Analysis", and I believe his textbook on hierarchical modeling: if a friend argues with you about this, I think those would be good places to turn to for help.
There are some who like to divide by twice $SD_x$, on the basis that this makes discrete explanatory variables and continuous explanatory variables directly comparable^[See Grueber *et al.* (2011; Journal of Evolutionary Biology 24, 699–711) for a discussion of this, particularly in the context of model averaging. The definitive paper is probably Gelman 2008; Statistics in Medicine, 27: 2865–2873.]. By all means, go nuts and do so (all my students do!), and you would (of course) be correct. I'm not necessarily sure I agree that it is philosophically possible to directly discrete factors and continuous factors (they just feel... different to me) but the math checks out so do what you wish.
So, from what we've learned, what is the next-best way to compare models after picking a sensible hypothesis and comparing that with a sensible null? Scale your variables, use model-averaging to get a set of standard effect sizes, and interpret those in terms of relative importance. Do this correctly, and you won't go wrong.
## Exercises
Today's data come from `boardgamegeek.com`^[Downloaded using code from <https://github.com/willpearse/boardgamegeek>—you can only imagine how exciting I am in my free time.]. They describe ratings given to board games by people on the Internet (`rating`), and information about each board game such as the `year` in which they were published, their `complexity` (how difficult they are to learn, rated out of 5), whether the game was funded through crowd-sourcing website `kickstarter`, and whether the game is co-operative (`co.op`; whether the game has a single winner at the end). As ever, below is some code to load today's dataset into `R`, critically and importantly there are two versions of this dataset: one is simplified and the other is more complex and has more variables. If you have attempted the extension/optional material, there are additional datasets that I also give instructions for the loading of below.
```{r}
# Load the simplified board game geek data in like this
data <- read.csv("bgg-simplified.csv", as.is=TRUE)
# Plot the data out like this - see if you can figure out what the ifelse is doing
with(data, plot(rating ~ year, col=ifelse(kickstarter, "red", "black"), pch=20))
with(data, plot(rating ~ complexity, col=ifelse(kickstarter, "red", "black"), pch=20))
with(data, boxplot(rating ~ kickstarter))
with(data, plot(rating ~ co.op))
# Load the complex board game geek data in like this
# (what does the row.names bit do, do you think?...)
data <- read.csv("boardgames-huge.csv", row.names=1, as.is=TRUE)
# Plot the data out like this
with(data, plot(rating ~ year))
# Here's how to fit a big model
top.model <- lm(rating ~ year + I(year^2) + max_players +
log.play.time + two.players + three.players + weight +
mech.pca.1 + mech.pca.2 + mech.pca.3 + kickstarter +
log.owned + log.n.ratings + co.op,
data=data, na.action="na.pass")
# ... what does the I(year^2) do? Why couldn't we write year^2?
# 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. Now you are ready to work with the simplified form of the data data.
a. Fit a regression to the data to see whether there is a 'cult of the new': whether newer games are rated more highly than older games.
b. Fit a multiple regression to the data to determine whether people more-highly rate games that were funded through kickstarter and are newer. Don't fit an interaction term in this model.
c. Test whether there is an interaction between age of game and whether it was funded through kickstarter.
d. Michael approaches you and says that board gamers are strange people who only like absurdly complicated games and care about winning more than anything. Using what you've already learnt about gamers' ratings of games, and the `co.op` and `complexity` variables in the dataset, fit a model to test their hypothesis. Optionally, but encouraged: flee from Michael, because surely you have learned from our first lecture to be wary around him.
3. Now you are ready to work with the more complex form of the board-game dataset.
a. Fit the standard regression above to these data. What do you think is going on?
b. Use model averaging techniques to figure out what's going on with this model (don't forget to use `sw`). Do you find it easier to interpret what's going on with these data now?
c. Scale the variables to contrast relative effect sizes. What is the most important factor in explaining a game's rating?
4. Today's code challenge is all about making it marginally easier for you to take coefficients from models calculated with standardized data and then, if you want, back-transform those coefficients. If you're not that great with algebra, take @eq-standardize and re-arrange that to get the right answer. A pen and paper will work wonders, and I promise you that if it seems awful you've likely done it wrong...
```{r}
#| eval: false
back.transform <- function(coefficient, orig.x){
sd.orig <- ____(x)
___.orig <- mean(x)
rescaled <- coefficient * ___ + ___
return(rescaled)
}
```
5. (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.
6. (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.
