# Hierarchical (Bayesian) models
In this session we're going to learn about hierarchical models, which are an important family of models that even crop up in machine learning. We can fit them in a frequentist framework, but only in special cases or through mathematical tricks such as so-called *mixed effects models*. Thus hierarchical models will form the perfect reason to introduce you to Bayesian frameworks. People often misinterpret me when I say that I hate mixed effects models: I don't hate them because they're complicated, I hate them because their complexity arises from their being statisticians' tricks to avoid soaking up degrees of freedom while estimating parameters^[See Pearse *et al.* (2014; Modern Phylogenetic Comparative Methods) for a review of a subset of them and Pearse *et al.* (2015; Bioinformatics) for software implementing them. So, to be complete, I don't hate them so much that I won't use them when they're appropriate, and they are sometimes appropriate (...but only when I say so...)]. I mention this because, ultimately, I think they're straightforward to understand if you know "the secret": they're all a big trick, and have no serious philosophical underpinnings. If you understand that, then it makes sense why some people have so much trouble defining what should(n't) be a random effect: there's no rule, because mixed effects models break the rules and are a trick. It also makes sense why it's so hard to get $p$-values for statistical significance from them: there's no 'official' way to do it, because they're a trick and so all the usual rules no longer apply. If you become interested in hierarchical modelling, I would recommend two excellent resources: the GLMM wiki (<http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#model-specification>) and Andrew Gelman's textbook on hierarchical modelling. Sadly, by the end of the Gelman book you won't be using mixed effects anymore, and will instead focus on Bayesian hierarchical models. Oh well.
## Mixed effects models
### An informal introduction to random effects
Imagine you are interested in the correlation between the amount of exercise a person does and their mood. You know that you can't study only one person, because that person might be an outlier and show an unusual response, so you measure multiple people. Just for fun, let's simulate some data that might correspond to this scenario:
```{r}
# Simulate some people with varying responses to exercise
person <- letters[1:20]
av.mood <- rnorm(length(person), mean=5)
exercise.effect <- rnorm(length(person), mean=1)
# Simulate some data from those people
mood <- as.numeric(t(matrix(exercise.effect,ncol=1) %*% matrix(0:1, ncol=2)))
mood <- mood + rep(av.mood, each=2) + rnorm(length(20*2))
# Reformat explanatory data to match
person <- rep(letters[1:20], each=2)
exercise <- rep(c(FALSE,TRUE), 20)
# Make a pretty plot where each person is a different color
cols <- setNames(rainbow(20), letters[1:20])
boxplot(mood ~ exercise, names=c("no exercise", "exercise"))
points(exercise+1, mood, col=cols[person], pch=20)
```
The code is a little complicated because it uses matrix multiplication (`%*%`), but the general idea is we sampled 20 people, measured their happiness after we made them exercise and after we forced them not to^[Please try not to think about how you would actually do this; just focus on the stats.], and then recorded their mood. We know that, on average, people are happier when they exercise (on average they are 1 mood point happier; see the third line), but we also know that people's average mood varies (see line two). Now we've got this data, can we detect an effect of exercise on people's mood?
```{r}
# Can we see an effect of mood?
summary(lm(mood ~ exercise))
#... ah, we can't. Maybe a better model would?
summary(lm(mood ~ exercise + person))
#...ah good, it can. But wait, how many parameters?
```
When I simulate it, we can't, but bear in mind that this is random data and so you might be able to when you run this example. We can, of course, when we include each person's overall mood as a factor (that turn this into a paired t-test, if it helps to think of it like that, but it doesn't do so literally) but there's a problem. When we estimate each person's mood, we estimate, overall, 21 parameters from only 40 datapoints. The general rule of thumb is that you need 10 datapoints for each parameter estimated; it works with this simulated data because I kept everything simple, but in more realistic simulations we still wouldn't be detecting an effect. Indeed, in more realistic simulations, we sometimes can't fit the model at all.
The proper way to solve this problem is to fit a *hierarchical model*. To do this we accept that each person's mood is, itself, drawn from a distribution (in this case a normal distribution with a mean of 5 and standard deviation of 1). We then model the parameters of the distribution from which everyone's average mood is drawn; we call these *hyper-parameters* because they sit 'above' everyone's moods and determine them. Indeed, we can use this to argue that while we're now estimating 22 parameters (20 mood parameters and the two hyper-parameters), we're not really soaking up 22 degrees of freedom—the hyper-parameters show us that the people's moods aren't independent, and they're to some extent explained away by the two hyper-parameters we've estimated. It turns out that the mathematics of all this are very, very complicated. So complicated, in fact, that it's not trivial to agree how many degrees of freedom we really are estimating^[A little more on this later in the chapter], and in only the simplest cases can we find an analytical solution for what the parameters in this model should be^[Loads on this next time].
So some bright spark came up with a solution: *cheat*. First, absolutely 100% promise that you will never, ever look at the estimated parameters for the average moods^[The only exception to this is that if you are only interested in estimating the moods themselves, and don't want to generalise from your model to any other data, you can look at the estimates. This is useful in a handful of cases where the mathematics behind such random effects allow for better estimation of estimates than fixed effects. Sadly, these rare cases are used to justify treating random effects as if they were fixed effects for convenience's sake, and it's *always wrong*. None of the techniques I am teaching you in this class make it justifiable to look at the random effect estimates: if you ever do it, you are wrong, and you must stop. It's unethical, and frankly stupid, to do otherwise.]. Second, assume that the effect of mood is essentially random—that the overall mean of the distribution from which the moods were drawn is 0. This isn't as bad an assumption as it sounds, because if we estimate an overall intercept/mean, all of the average moods are really just offsets, or contrasts, from that overall mean. Finally, if we do all that, then the only parameter we have to worry about is the standard deviation of the moods. So estimate that standard deviation and you're done: we can account for the effect of individual-level variation while only estimating one parameter.
Thus we have modelled `person` as a *random effect*, and can continue to model `exercise` as a regular effect, which in this context is called a *fixed effect*. Mixed effects models are so-called because they mix fixed and random effects. Whole blog posts, and even articles, have been written about whether you should treat a variable as a random or fixed effect. It boils down to this: if you want to look at the values of the parameters you are estimating and use them for something, then you need to fit that variable as a fixed effect. If you want to generalise something about that variable to another system, or make a prediction using it, then you're essentially looking at that value and so, again, you need to fit it as a fixed effect. Otherwise, you essentially don't care about that variable but you need to account for its effect: thus you can model it as a random effect.
People often describe random effects as a good way to handle *pseudo-replication*. Pseudo-replication is a fancy word for saying that a biasing factor makes it look as though you've sampled something more often than you have. For example, if I said I had 100 samples of mood, but they were all taken from the same person, you'd say my 100 samples aren't totally independent: they are pseudo-replicated. Adding a random effect that matches onto the level of pseudo-replication is one way of accounting for that.
### Hands on with mixed effects
```{r}
library(lme4)
mixed.model <- lmer(mood ~ exercise + (1|person))
summary(mixed.model)
# Once you run the next line, the output of the line above will change
library(lmerTest)
summary(mixed.model)
```
The code above fits a mixed effects model to the data we've simulated, and then displays its output using two different packages (hence the comment). The model formulation `(1|person)` means "estimate the intercept (parameter 1) and then, as a random effect varying about that, fit the variable person". A nice thing about mixed effects models is you get the random effects' variance ($\sigma^2$; its square root is standard deviation $\sigma$). It gives you an indication of how much variation there is among simulated people; it also makes it possible to *partition* variance among hierarchically-nested groups. Something like this:
```{r}
# Simulate a weird, useless higher-level group
higher.group <- ifelse(person %in% letters[1:10], "one", "two")
# Fit new model
silly.model <- lmer(mood ~ exercise + (1|higher.group/person))
# Perform variance composition analysis
vars <- as.data.frame(VarCorr(silly.model, comp="variance"))$vcov
names(vars) <- c(names(VarCorr(silly.model)),"residual")
(vars / sum(vars)) * 100
```
We could make a statement now along the lines of "roughly 85% of the variation in mood is attributable to the person who was sampled, and essentially all of the rest is residual variation", although your precise values will vary depending on how your data ended up being simulated.
You may remember that, when we simulated our data, we allowed each person's response to exercise to vary. We can, if we wish, account for this in a mixed effects model too, although `lme4` will complain if we don't have a greater replication of data (even mixed effects models have limits!). So, let's simulate some more data and fit a more complex model:
```{r}
# Repeat simulation with more data
mood <- as.numeric(t(matrix(exercise.effect,ncol=1) %*% matrix(0:9, ncol=10)))
mood <- mood + rep(av.mood, each=10) + rnorm(length(mood))
# Reformat explanatory data to match
person <- rep(letters[1:20], each=10)
exercise <- rep(0:9, 20)
# Fit new model
complex.mixed.model <- lmer(mood ~ exercise + (exercise|person))
```
This is quite a complex mixed model, but it's also of the form that most people want to fit and so I think it's useful to present it to you. The random effect term `exercise|person` means "allow the intercept and the slope to vary across each person", and means this is something called a "random slope" model. Please ignore people on the Internet who state that such models should be specified as `1+exercise|person`; the two formulations have always been identical and there's no need to make the syntax more complex than it already is. What makes this model so complex for the computer is that `exercise` appears as a fixed factor as well: this fixed factor represents the overall impact of exercise on mood. Again, despite what you will read on the Internet, it's totally valid (indeed, quite normal) for terms to appear in both the fixed and random parts of a mixed effects model.
### Significance of fixed and random effects
I'm afraid it's impossible to talk about testing the statistical significance (*i.e.*, getting a $p$-value) of a mixed effects model without discussing how mixed effects models are numerically fit. I can't appeal to all the statistical and mathematical niceties that I employed earlier in the course because we left all of those behind when we applied the mathematician's trick. Sorry; strap in.
When fitting a mixed effects model via maximum likelihood, there are two approaches: you can estimate your fixed effects first and then *then* your random effects, which is called *Maximum Likelihood* (ML) fitting in this context. Or you can estimate the random effects first, and then estimate your fixed effects, which is called *Restricted Maximum Likelihood*, or sometimes *Residual Maximum Likelihood*, and is abbreviated to REML. If you want to say something about the random effects, you've got to fit under ML, and if you want to say something about the random effects, you've got to fit under REML. Each is calculated in the context of the other, *and you cannot estimate the significance of both at the same time*.
Another word about pseudo-replication. Mixed effects models are a very good way of accounting for unequal or repeated sampling. Let's say that you asked the same person about their mood multiple times: clearly that's a form of pseudo-replication, because each mood estimate is not independent. Before a mixed-effects modelling course, you might have averaged across those mood estimates in order to remove that pseudo-replication, but you would also be throwing away lots of information. A hierarchical model—of which mixed effects models are the only example we know of right now—let you measure and account for that pseudo-replication using random effect terms. So if someone every tells you to 'average out' something, take a minute and think whether it might be best to fit a mixed-effects model or something of that ilk.
#### Random effects
I rarely see the point in estimating the significance of a random effect. If you're using your random effects to handle pseudo-replication^[...and, honestly, if you're not then why aren't you fitting your variables as fixed effects?] then I can't see a defensible argument as to why you should have to prove that your data are pseudo-replicated in order to account for it. Philosophically, they either are or they aren't: choose, and live with it. The only exception to this that I can see if where you are interested in the variance component of the random effect as a variable of interest in-of-itself. An obvious example (to me) is Phylogenetic Generalised Linear Mixed Models, where we use random effects to represent the hypothesised rate of evolution of a trait across a phylogenetic tree. Brownian motion evolution has a single parameter, $\sigma^2$, which can be viewed as the variance of a multivariate normal distribution from which we draw intercepts/slopes for species. We care about that variance: if it's indistinguishable from zero, then we don't think Brownian motion evolution is going on. So I think it makes sense, in that context, to test whether that parameter is zero^[I might add that I don't like it, and am working on better alternatives, but this is all we have right now.]. Such models are also relatively common in population genetics and genomics, and the so-called 'animal model' is a prime example.
All this to say that, apart from very exceptional cases where quite a lot of argument and mathematical hand-wringing is required, I don't see why you should be testing the significance of random effect terms. If you do find yourself doing something that, to be clear, I really don't think you should be doing, your best bet is a *Likelihood Ratio Test*, which is a generalisation of the $F$-tests that I've already taught you. In other words, use our old friend the `anova` function and you should be fine:
```{r}
random.slopes <- lmer(mood ~ exercise + (exercise|person), REML=FALSE)
random.intercepts <- lmer(mood ~ exercise + (1|person), REML=FALSE)
fixed.only <- glm(mood ~ exercise)
anova(random.slopes, random.intercepts)
anova(random.intercepts, fixed.only)
```
Notice that I re-fitted the models under REML; some packages will warn you if you forget but not all. Notice also that, we have to specify that we need to perform a $\chi^2$ test, which, if you did the GLM extension exercise you will recall vividly, and can just about see if you cast your mind waaaaay back to lecture 1, are really essentially the same as $F$-tests. In these cases the more complex random effect structures explain more variance than the simpler/absent ones, and so we would conclude that they were 'significant'. I air-quote significant because the variance term is bounded to be positive^[It's a variance; you can't have negative variances.] and all ML-based tests perform a bit poorly when faced with such bounded distributions of parameter values.
I know some people who fit things that aren't significant when treated fixed effects as random effects in order to 'make them' significant. These people are unethical: I've already gone through the problems with assuming something is a random effect. Just don't do it. People tend to do things like this in order to publish Nature papers; do consider that you may well get found out and it's not a good look.
#### Fixed effects
The short answer to "how do I test the significance of a fixed effect" is that you can't. I appreciate that you won't be very satisfied with that answer, but it's the truth. You can't. Sorry. Entire courses are structured around mixed effects models (we're only here to get to the next page, really), but that is the truth. If you don't believe me, I encourage you to read Doug Bates' (in)famous discussion of this on the `R` mailing list (<https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html>). If you still don't believe me, load up the help file for `lme4` and check out the help file entry for `pvalues`. It is damning that the world experts in mixed effects models (Doug Bates, Ben Bolker, and friends) literally refuse to write a $p$-value function into the package that everyone else relies on. The 'standard' way to get $p$-values for fixed effects is to load a package written by someone else (they're listed in that help file entry, by the way) and use their 'approximation'. The alternative is to use something called 'MCMC sampling', which is essentially what I'll teach you next time.
The programmers at SAS created an approximation to an $F$-test that someone ported into `R`; you can get to it like this:
```{r}
mixed.model <- lmer(mood ~ exercise + (1|person))
summary(mixed.model)
# Once you run the next line, the output of the line above will change
library(lmerTest)
summary(mixed.model)
```
It's quite an amazing piece of code, really: it masks `lme4`'s functions to give you identical output save the addition of a $p$-value. Ultimately, of course, it's a nonsense: the fundamental problem with random effects is we are only *pretending* that they soak up only one degree of freedom to estimate their variance. Thus we don't know how to calculate an $F$-statistic, because we don't really know how many degrees of freedom we're dealing with. There can be no 'correction'; we literally have no idea, philosophically, what to do—we're just making things up. Ironically, calculating a random effect's significance is easier because we're just looking at the variance term. That we know how to deal with; everything else, we don't.
It is possible to perform AIC-based model averaging of mixed effects models, and so you can avoid the $p$-value problem that way. Of course, AIC involves a term for number of parameters estimated, so the whole thing remains something of a farce, but if you keep the random effects constant it's at least an informative farce. Standardising your variables will make their relative effect sizes comparable, and is a trick I've pulled on several occasions. I'll finish as the `lme4` help file for $p$-values does, and direct you towards PhDComics: <http://phdcomics.com/comics/archive.php?comicid=905>.
### (Extension) Model criticism of mixed effects models
(Before reading this section, ensure you have read the extension section on model criticism.)
It is impossible to define an $r^2$ for mixed effects models because we don't really know how many degrees of freedom they soak up^[The same is not true of a full hierarchical model, where it is possible to have a full range of non-integer degrees of freedom to account for this non-independence. If you're interested, read Gelman & Pardoe's "*Bayesian measures of explained variance and pooling in multilevel (hierarchical) models*" (2006; Technometrics).]. A definition does, however, exist: a *marginal* $r^2$ that measures the effect of the fixed effects, and a *conditional* $r^2$ that measures the effects of both the fixed and the random effects (see Nakagawa & Schielzeth, 2012, Methods in Ecology & Evolution where they were defined). You must never interpret one of these metrics without the other, but, depressingly, that doesn't stop people from doing so. Bluntly, I think it's all a lot of codswallop: if you look at the definitions of these metrics, they're quite patently *not* $r^2$ values and defining something that is almost-but-not-quite an $r^2$ and then calling it an $r^2$ does not, in fact, make it an $r^2$. If they were $r^2$ values, then you wouldn't need two of them! Regardless, they are useful to know about, but make sure you treat them with a judicious pinch of salt.
Always plot both your raw data against your predictions to see if your model is behaving sensibly. Such diagnostic plots are a special case of something called *posterior predictive checks*, and are informally called "making sure your model doesn't predict nonsense". Never believe anyone who says their mixed effects model is too complex to plot the raw data: that person either doesn't know what they're doing, or does and is trying to hide something. People will often talk about checking 'within group variance'; what this means is plotting their data out across groups against the predictions from their model to see if it's capturing the variation within and among groups. This is perhaps the only time you're allowed to look at the estimated random effects (using something like `ranef(model)`). If the predictions from the random effects match nicely onto the raw data, then your model has accounted for the variation in the random effects well. If it doesn't, it hasn't, and you should fit a more complex random effects structure. Given that you should be fitting as complex a random effects structure as you can, for the philosophical reasons I outline above, please don't use these plots as an excuse to 'drop' random effect terms. If you suspect pseudo-replication, the suspicion is enough, and getting rid of something to make your fixed terms significant is very bad practice indeed!
## Fully hierarchical (Bayesian) models
Now it's time to overcome those problems I've been banging on about by moving to a fully-specified hierarchical model using a Bayesian approach^[To be clear: mixed effects models are hierarchical models too. What I mean by 'fully specified' is that we will no longer be applying cheap tricks that permit estimation at the cost of disallowing interpretation of those estimates.]. We're going to be covering a lot of ground so it's important that you keep two things in mind: the definition of a Bayesian probability that we encountered in chapter one, and what we learned about mixed models in the last section. Now we're going to be removing the complications of frequentist mixed models and 'going full Bayesian': estimating the importance of different terms in our model with full confidence. An important thing to bear in mind about Bayesian models is that, really, they're not as complex as everyone makes out. The Bayes formula is not that complex (although, as we'll see, there are some common misconceptions about it), although Markov Chain Monte Carlo (MCMC) can seem a little weird. My advice is to not worry about where the mathematical proofs of MCMC come from: you've accepted, quite happily, lots of things about frequentist statistics without formal proof^[Regression equations go through the mean, that $F$-statistics can be standardised according to degrees of freedom, that the maximum likelihood estimate is the same thing as the estimate that minimises the sum of the squared errors... If I stop now do you promise not to be confused by what you already know?]—my advice to you would be to let MCMC be the same and "just trust me".
### Practical Bayesian inference
Armed with our new-found experience, let's formally define the concept of a **frequentist** probability:
$$
p(data|model) = p(x < a \cup x > b) = \int_b^a \text{pdf}(x) dx
$$ {#eq-likelihood}
Where $model$ is some statistical model, $|$ means "given", as in "given we're dealing with some $data$", $a$ and $b$ are possible parameter estimates, $\cup$ means "and"^[Technically this is a 'union', but the meaning is the same.], as in "both of these things are true", $\int$ is an integral sign and means "area under a curve", *pdf* is a *probability density function* (@eq-normal is the probability density function of the normal distribution), $dx$ means our probability density function (and so area under our curve) is calculated across $x$, and so $x$ represents all the possible values of our coefficient. Phew! A lot of words, but, ultimately, it's nice to have everything written out at last, wouldn't you agree^[Oh, alright, one more term if you want it: a formal definition of the "likelihood" (*i.e.*, the thing we tend to maximise) would be $\text{pdf}(\theta)$.]? In words, this means exactly what we're already familiar with: the probability that we would some data like our draws from a particular distribution (*pdf*) is equal to the area under that distribution's curve ($\int$) between two points on that curve ($a$ and $b$). When we were using this to measure how surprised we were by a particular $t$ or $F$ value, we set $a$ and $b$ to be something like $1.96$ and $\infty$ (infinity): we were interested in the right- and/or left-hand tails of the distributions.
I write all of this out to highlight two things. First, frequentist probabilities can, and indeed often do, involve conditionals (the $|$): Bayesians don't have a monopoly on it^[The Internet is filled with "Bayesian solutions to the Monty Hall problem"; there are such solutions, but the solution that's often given is actually frequentist. Updating your probability/likelihood in the face of new information (*i.e.*, data) is a frequentist thing too; if you find such an example online and become confused, come ask me about it.]. Second, the formal **Bayesian** definition of probability is actually much simpler, in many ways:
$$
p(model|data) = \int_b^a posterior(x) dx
$$ {#eq-bayes_prob}
Two new things here: we're calculating the probability of our *model* given the data we have in hand, and we're working with something called the *posterior distribution* now. This is much cooler, actually, because as well as being simpler our probability is actually the thing we wanted to know all along: how likely we are to have got the right model/coefficients. This means we can say things like "it's ten times more likely that our two treatment groups differ than they don't" (*i.e.*, our probability of a difference is $\frac{10}{11} \approx 0.91$), which we can't say in frequentist land. The only snag is we need to define the posterior distribution:
$$
p(model|data) = \frac{p(data|model) \cdot p(model)}{p(data)} = \frac{likelihood \cdot prior}{normalising\ constant}
$$ {#eq-bayes_theorem}
@eq-bayes_theorem is essentially the *Bayes Rule* and I've given it in both words and equations above. Notice that it makes use of @eq-likelihood, as well as a *prior belief* about both our model formulation (that our data are from the normal distribution, for example) and parameter estimates (that our normal distribution probably has a mean of 5, but definitely doesn't have a mean of $-1000$, for example), and a normalising constant that's essentially the probability of our data. We'll cover each of these concepts in their own sections below, but first let's just think about what a posterior distribution represents. Instead of having an *estimate* with *uncertainty*, as we have in maximum likelihood, you now have a distribution whose centre is essentially your estimate and whose variation is the uncertainty. What's magical about this distribution, once again, is taking an area under it corresponds to calculating the probability that your estimate is correct. Thus you might have a posterior mean/median of 5, and 95% of the *posterior density* between 4 and 6. This means your mean estimate is 5, and you're 95% certain that your estimate is between 4 and 6^[Such intervals are called *credibility intervals*, and are often abbreviated to CI, in order to reassure frequentists who are desperate for '95% CIs'. Of course, these CIs are much more meaningful and intuitive!]. Isn't that great?
### A brief digression into priors
The concept of a 'prior belief' is something that gets people really upset, for reasons that aren't always totally logical. Bayesian statistics requires the specification of a set of prior beliefs about what parameters values might be reasonable for your model: you literally specify relative probabilities for each possible parameter value in your model in the absence of data. Let's say, for example, you were trying to estimate the variance of a normal distribution ($\sigma^2$): a sensible prior would be to say that, in the absence of any information, you think that parameter could be, with equal probability, somewhere between $0$ and $\inf$. This required two decisions: the specification of a distribution—in this case the Uniform distribution—and parameters for that distribution. The specification of the distribution should be determined by the kind of parameter you're using: each distribution's parameter has a *conjugate prior distribution* that you should use to ensure your inference isn't biased^[...there's actually not always a conjugate prior, and it's quite often not the distribution you would think, but I don't have time to explain all the details.]. Our variance parameter here is a good example: you can't have a negative variance, so we need to pick a prior that doesn't allow negative values. The next question is what parameters to use for that prior distribution, and there is no hard-and-fast rule for how to do this other than it should reflect your own prior beliefs. If you don't know, and I would argue even if you do, then make your prior nice and wide (as in this example, which is infinitely wide), and then you're fine^[...it can also be more complex than this. In cases where your model doesn't work well, you may need to tighten your prior to essentially 'nudge' your model in the right direction. But see my advice in the text.].
People tend to get very concerned about priors because there truly is no answer to the question "what priors should I use?". So bear the following two things in mind. (1) There are two things that determine your posterior distribution (your model output): the data and the prior. If you change your prior and your posterior/answer changes, then that suggests your prior is "swamping" your data, and that you can't be certain given your data. Thus, if you're worried, change the prior drastically^[You must, of course, make sure that your prior doesn't preclude what you found before by design: if the 'true' mean is 7 and you set a prior with a 0 probability for 7 then, of course, you won't get the right answer.] and see what happens: if nothing changes, you're fine^[...unless you can make an argument as to why your prior is correct. If you can make such an argument, good for you, but I can rarely convince myself of such stuff.]. (2) Remember that likelihood are usually expressed in log-space and the same is true of priors. Thus, when calculating a posterior distribution, you add up as many log-likelihood estimates as there are datapoints and only one number for the prior (regardless of the number of datapoints). Thus the prior is a tiny, essentially inconsequential, part of your model.
### Normalising constant
Without getting into the nitty-gritty of it, the normalising constant is the reason that everyone left Bayesian statistics alone for so very long. It turns out it's difficult, if not essentially impossible, to calculate the probability of your data. Indeed, it's not just numerically difficult (*i.e.*, would take a lot of computing time), it's also analytically difficult—in many cases, we're not sure if we'll ever be able to figure out the true mathematical solution.
Luckily, there is a way forward, using *Markov Chain Monte Carlo* (MCMC). It turns out that by simulating (that's the *Monte Carlo* bit, named after the city with all the casinos as it's a random process) one or more *Markov Chains* (a process where the only thing that determine your future value is your last value) for a long time, we will eventually *converge* on the true posterior distribution without ever having to solve the complicated maths. There are many kinds of *MCMC algorithm*, but here is an outline of one of the simpler ones—the *Metropolis-Hastings algorithm*:
1. **Setup.** Pick random starting values for each of your model parameters—normally by drawing from your prior distribution.
2. **Proposal.** Draw, from some random distribution called the *proposal function*, new parameter value(s).
3. **Calculate.** Calculate the *acceptance ratio* $\alpha$: divide the likelihood of the model with the proposed parameter values by the likelihood under the old values.
4. **Accept or reject.** Accept the new parameter values, at random, with probability $\alpha$.
5. **Iterate.** Go back to step 2 (proposal).
This process repeats until the end of time or you tell the computer to stop. Given some assumptions about your proposal function (which are remarkably lenient), in principle this process should always *converge* on the *posterior distribution*—you will end up sampling a series of parameter values whose histogram is indistinguishable from that of the posterior distribution.
The calculus as to why this works may seem like magic, and perhaps it truly is. Regardless, there are lots of different algorithms that should, in principle, converge upon the posterior distribution. The one we will be using today is called *Hamiltonian Monte Carlo*, and is essentially the same as described above with the addition that the Markov Chain (the process of picking new states with each iteration) has a degree of *momentum*. This is mathematically identical to an ice-skater moving across ice, where the ice is the chain's position in parameter space and relative height of the ice is the relative likelihood of different positions in space. As the chain 'skates' downhill from regions of bad likelihood (which are *higher* regions of ice) it picks up momentum, allowing it to find the better regions of the ice-rink more quickly. If this seems like nonsense to you, don't worry, and just focus on understanding the general process above^[If you are programmatically inclined, you might want to look at a toy MCMC implementation I wrote some time ago: <https://github.com/willpearse/austin>. Go back through the commit log; it starts with a basic MCMC and works up to hot/cold Multiple Chain MCMC.].
As with everything in stats, there is no way to tell whether your MCMC has converged onto the posterior distribution. The clues we use to figure it out are:
* Running multiple chains and checking to see if they have *mixed*; *i.e.*, their distributions are indistinguishable.
* Plotting the posterior probability and/or likelihood over the course of the MCMC run. It should shoot upwards and then randomly dart around ("the hairy caterpillar"), indicating that you've hit the posterior distribution. Similar plots of coefficient values 'converging' on some sort of reasonable value can be made.
* The posterior should, ideally, not look exactly like the prior distribution. If it does, make sure you haven't mis-specified your model somehow. If you are using `BUGS`, `JAGS`, or `STAN` (we are not today—we are using something derived from them) then check your code as you've probably got a bug.
It is absolutely vital that you should *never* make inferences from a posterior distribution that either hasn't converged, or contains portions of the chain that haven't converged. MCMC is perfectly mathematically valid, but the whole house of cards falls down if it hasn't converged. The section of the chain before convergence is called the *burn-in*, and should be removed before any formal inference. While there are metrics that tell you whether or not it looks like a chain has converged, those bring us right back to the whole Gödel's incompleteness thing: they have assumptions, and they need to be tested somehow. The bottom line is: it's never a "near miss", and so you can always tell easily if a chain hasn't converged. You should always cut off more as burn-in than you think you need, and 10% of your chain is a good minimum (processing time is cheap; run the chain for longer to get there if needs be).
### Hierarchical models formalised
In theory, all that remains now is to formalise a hierarchical model a little more formally. We could, in principle, charge ahead without this but I think it's nice to be thorough.
Imagine we have a continuous response variable ($y$) that we're modelling as a function of the interaction between continuous and discrete variables:
$$
\mu = a_i + b_i.x
$$
$$
y \sim Normal(\mu, \sigma^2)
$$
Where $y$ is our response variable, $\mu$ is our predicted response, $x$ our continuous explanatory variable, $a_i$ is the 'intercept' estimate for each level ($i$) of our discrete variable (*i.e.*, the mean of the response variable for each level), $b_i$ is the 'slope' estimate for each level ($i$) of our discrete variable, and $\sigma^2$ is the variance (the amount of error, essentially) we would expect. The only thing that's new here is the way it's written: I've used subsetting ($_i$) to indicate estimates for each level of our discrete explanatory variable.
To fit this kind of model in regular `R`, we'd write something like `lm(y ~ x * discrete)`. However, we could make this a *hierarchical* model by defining the following:
$$
a_i \sim Normal(\alpha, \sigma_\alpha^2)
$$
$$
b_i \sim Normal(\beta, \sigma_\beta^2)
$$
Now we've defined four *hyper-parameters*: $\alpha$, $\beta$, $\sigma_\alpha^2$, and $\sigma_\beta^2$. These parameters define how we think the slopes and intercepts in the definition of $\mu$ are, themselves, defined. They're called *hyper*-parameters because they sit 'above' the parameters in our model; they're the parameters of the parameters. In a Bayesian setting we would have to define priors on these hyper-parameters, which we could call hyper-priors if we wanted to confuse people. Notice that I didn't have to define two sets of hyper-parameters; I could have only had the $\alpha$ or only the $\beta$ terms in there and it would still have been a valid hierarchical model. Notice that it's quite common to use subscripts and/or Greek letters to distinguish between hyper- and regular parameters.
A hierarchical model is defined as a model with hyper-parameters in it and, just like we saw in mixed effects models (which are hierarchical models), you can nest the equations to have hyper-hyper-parameters (just call them hyper-parameters, though). We can define mixed effects models in this framework; the difference is that we have to explicitly introduce terms for the 'overall' intercept and 'overall' slope. For example, a mixed effects version of the above would be:
$$
\mu = a_i + a_0 + (b_i+b_0).x
$$
$$
y \sim Normal(\mu, \sigma^2)
$$
$$
a_i \sim Normal(0, \sigma_\alpha^2)
$$
$$
b_i \sim Normal(0, \sigma_\beta^2)
$$
Where $a_0$ and $b_0$ are the overall/mean intercept and slope terms, respectively, and we have now removed the $\alpha$ and $\beta$ terms. If mixed effects models seem like a bit of a trick, in the sense that we've just moved a hyper-parameter 'down' to a regular parameter, that's because they are...
### Hands on with `rstan(arm)`
The good news is that fitting Bayesian hierarchical models is a lot easier than it used to be. You need to install a package called `rstanarm`, which can take a very long time to install, but once it's loaded you'll all ready to go. It will even choose your priors automatically for you, which is fantastic. Let's start off by simulating the same data as we were using last time, and then fit a Bayesian hierarchical model to the data.
```{r}
# Repeat the last simulation from the last session
person <- letters[1:20]
av.mood <- rnorm(length(person), mean=5)
exercise.effect <- rnorm(length(person), mean=1)
mood <- as.numeric(t(matrix(exercise.effect,ncol=1) %*% matrix(0:9, ncol=10)))
mood <- mood + rep(av.mood, each=10) + rnorm(length(mood))
# Reformat explanatory data to match
person <- rep(letters[1:20], each=10)
exercise <- rep(0:9, 20)
# Load package
library(rstanarm)
# Fit model
model <- stan_glmer(mood ~ exercise + (1|person))
```
You'll see that it takes a fair bit longer to run the model and its MCMC chain, but the syntax is identical to that of `lmer`. Notice that it runs four separate chains by default (each of length 2000), and throws away the first 50% of each chain as a burn-in. You can play around with these settings if you wish, of course.
```{r}
# A shorter set of coefficients
model
# More comprehensive and useful coefficients
summary(model, probs=c(.025, .5, .975))
# "Caterpillar" coefficient plots
plot(model)
# MCMC diagnostics (see note in text)
plot(model, "trace")
```
The first kind of output is almost maddeningly terse, but recognisable from our mixed effects model days. You can see the median values of the parameters calculated across the posterior distributions: these are our 'best' estimates, and the *Median Absolute Deviation* (MAD) plotted for each parameter. You can interpret these exactly as we did last time for mixed effects models, right down to the error terms (the $\sigma$ terms from the random effects last time).
There are still no $p$-values in this output, but there are in the second—the are just hiding. By asking for the estimates of the 2.5th and 97.5th quantiles (and the medians) of each coefficient in our posterior distribution, we can see the range of coefficient estimates between which we are 95% sure our true value lies. In mine, for example, the slope of `exercise` has 2.5th of 0.6 and a 97.5th of 1.1. Thus I can say there is a 95% probability that the slope of exercise is between 0.6 and 1.1. I don't have to say I'm more surprised than I would be 5% of the time, or any of that nonsense: there's a 95% chance that the slope is somewhere between 0.6 and 1.1. How refreshingly simple. Even more excitingly, because random effects don't really exist in true hierarchical models^[I would argue they never did anyway, but moving on...], we are given estimates of the mood of each person and we are now allowed to do whatever we want with them. Indeed, those estimates have ranges of uncertainty too (the technical term is *posterior credibility intervals*, if you must know). Wonderful^[If this doesn't seem wonderful to you, please consider that I *wrote my own C++ MCMC sampler because I found all the alternatives too difficult*. This has become almost terrifyingly straightforward in the last few years!]!
The last two lines give us a plot of the coefficients and their credibility intervals, and then plots from the MCMC chain. For some reason `rstanarm` doesn't keep the burn-in phase of the MCMC chain in memory, so if you can see any kind of trend to these plots there's a huge problem. More details are given in the help-file entry for `plot.stanreg`. For all other model criticism problems, follow the guidelines in the mixed effects model section.
## Exercises
As ever, below is some code to load today's dataset into `R`. You're going to start off by working with it in 'mixed effects' mode, and then move onto analysing it using Bayesian hierarchical models. We then have a more theoretical exercise to test your understanding of Bayesian approaches, and then some extension simulation exercises for those of you who are interested in power analysis (estimating the power—$\beta$—of statistics in real-world settings).
```{r}
# Load the data in like this
data <- read.csv("climate-change.csv", as.is=TRUE)
# Plot the data out
state.cols <- setNames(rainbow(length(unique(data$State_abb))), unique(data$State_abb))
with(data, plot(human ~ temp.change, col=state.cols[State_abb], pch=20))
```
1. Your first task today is to go through the code above and make sure you can fit all the models described, and that the lectures make sense.
2. Today's dataset is an amalgam of many sources: a study on climate change perceptions (Howe *et al.* 2015; Nature Climate Change, 5(6), 596), the 2016 US Presidential election results (from Wikipedia), and NOAA's index of US temperature change (1901–2015; taken from the EPA website). Your task is to model how the percentage of surveyed people who think that climate change is happening (`human`) varies as a function of the observed temperature change to date (`temp.change`), recent voting records (`trump.percent`), and variation with and among states (`State_abb`). Note that there is an estimate in this dataset for (almost every) count within the US. I emphasise, strongly and clearly, that I am not ascribing any positive or negative connotations to anyone's voting intentions, I'm just trying to make a dataset that's slightly more interesting for you to work with than crabs that I have caught.
a. Fit a standard linear model (`lm`) to these data, modelling whether people think climate change is happening as a function of recent temperature change and the percentage of the state that voted for Trump. What do you find? Given the sampling of the data, do you think your model is reliable?
b. Fit a mixed effects model to your data, using a random effect to account for the repeated measurements within each state. How does this change your results?
c. Produce a plot showing your results. Interpret the graph, and explain whether you think the model is a good fit to your data.
d. Fit a Bayesian hierarchical model to your data, allowing each state to have its own intercept. What can you conclude from your model?
e. If you are feeling brave, fit a hierarchical model with variation in slope across states. How reliable do you think this model is?
3. Below is some code that is a toy implementation of MCMC to study the blood pressure data Michael collected for lecture 1's exercise. Test your understanding by filling in the answers below and carrying out the exercises.
a. I have highlighted the five steps of MCMC in the code comments below. Label those steps.
b. Generate a plot of the posterior probability across the iterations, and write down where the burn-in takes places. Note that the output from the MCMC is stored in a variable called `output`, and from following the plotting instructions at the end of the code you should be able to determine what each of the three columns in the output are.
c. Write out the prior specification and the form of the model that is being fit to these data (*i.e.*, find the response variable, the distribution is it being drawn from, and the parameters of that distribution).
d. What is the Bayesian estimate of the mean of these data? On the basis of this output, what is the probability that the mean of the data lie between 115 and 125?
e. Modify the priors to ensure that your prior specification isn't unduly affecting your answer.
```{r}
data <- c(120, 130, 110, 105, 121)
n.iter <- 10000
# The following two lines define the priors
mu.prior <- function(x) dnorm(x, mean=0, sd=10, log=TRUE)
sd.prior <- function(x) dnorm(x, mean=0, sd=10, log=TRUE)
# The following line defines the posterior distribution
posterior <- function(x, mu, sigma) sum(c(
dnorm(x, mu, abs(sigma), log=TRUE),
mu.prior(x), sd.prior(x)
))
proposal <- function(x) x + rnorm(1, 0, sd=.25)
# What step are the following two lines?
mu <- 90
sigma <-5
output <- matrix(NA, nrow=n.iter+1, ncol=3)
colnames(output) <- c("mu", "sigma", "posterior")
output[1,] <- c(mu, sigma, posterior(data, mu, sigma))
# What step is the next line?
for(i in 1:n.iter){
# What step is the next two lines?
new.mu <- proposal(mu)
new.sigma <- proposal(sigma)
post <- posterior(data, mu, sigma)
new.post <- posterior(data, new.mu, new.sigma)
# What step is the next line?
if(runif(1) <= 10^(new.post-post)){
mu <- new.mu
sigma <- new.sigma
post <- new.post
}
output[i+1,] <- c(mu, sigma, post)
}
# Plot out a column of the data like this
hist(output[,1])
plot(output[,2])
```
4. (Extension) Today's `R` challenge is an informal introduction to something called *power analysis*. You're going to simulate some data drawn from two different distributions, and see how varying the number of samples from each distribution affects your chances of detecting a difference between the two. This isn't something you would do with a mixed effects model necessarily, but it is something you might want to think about doing before you conduct an experiment. Fill in the blanks in the code below to make the simulation work.
```{r}
#| eval: false
# Setup simulation (look at expand.grid line as it's a neat trick)
difference.in.means <- 0.5
simulation <- data.frame(expand.grid(
number.samples=seq(5,100,by=2), number.replicates=1:10, p.value=NA
))
# Start simulation (with a 'for' loop)
for(i in 1:nrow(simulation)){
first.sample <- rnorm(simulation$number.samples[i], mean=0)
second.sample <- rnorm(simulation$number.samples[i], mean=____)
result <- ____(first.sample, second.sample)
____$p.value[i] <- result$p.value
}
# Plot the results out in an unhelpful way
with(simulation,
plot(p.value ~ number.samples, col=ifelse(p.value<0.05, "red", "black"), pch=20)
)
# Summarize the results a little
prop.signif <- with(simulation, tapply(p.value, number.samples, function(x) sum(x<0.05)/length(x)))
n.samples <- as.numeric(names(prop.signif))
# Plot results. Fill in your expectation of the power (from lecture 1) below
plot(prop.signif ~ n.samples)
abline(h=___, lwd=3, col="red")
```
5. (Extension) Now let's make the power analysis above even more complex. Simulations such as these are common for statisticians when faced with questions like "how many samples should I take given I know nothing about the system?". The answer is *always* "do a power analysis". "*To call in the statistician after the experiment is done may be no more than asking him to perform a post-modem examination: he may be able to say what the experiment died of*" (RA Fisher, 1938).
```{r}
#| eval: false
# Setup simulation (look at expand.grid line as it's a neat trick)
simulation <- data.frame(expand.grid(
number.samples=seq(5,100,2), effect.size=seq(.1,2,.2),
number.replicates=1:10, p.value=NA
))
# Start simulation (with a 'for' loop)
for(i in 1:nrow(simulation)){
first.sample <- rnorm(simulation$number.samples[i], mean=0)
second.sample <- rnorm(simulation$number.samples[i], mean=____)
result <- ____(____, ____)
____$____[i] <- result$____
}
# Summarize the results a little
prop.signif <- with(simulation,
tapply(p.value, list(number.samples,effect.size), function(x) sum(x<0.05)/length(x))
)
# Make some plots
n.samps <- unique(simulation$number.samples)
eff.sizes <- unique(simulation$effect.size)
image(prop.signif, x=n.samps, y=eff.sizes, xlab="____", ylab="____")
contour(prop.signif, x=____, y=____, add=TRUE, nlevels=5)
```
