1 Philosophy of statistics

In this first chapter, I’m going to try and give you an overview of the core concepts that define statistics. The only way this can possibly work is if you bear with me and promise to try and forget everything you think you have been told about statistics¹. It’s easy to get the wrong end of the stick with these concepts, particularly if you first encountered them some time ago, and often the easiest way to learn is to start again. If you’re coming from a discipline other than statistics, please also try and bear in mind that statistics is likely a larger discipline than your own subfield. This means that, in order to develop a serious statistical intuition and understanding, you’ve got to learn some fundamental vocabulary and concepts. Don’t approach this class with the mindset of trying to figure out how it will help you analyse the dataset you have in front of you right now: treat it as a grounding in the fundamental concepts that will allow you to learn anything you want about any dataset. If you’re impatient and just want answers, then fair enough, but please just bear with me…

1.1 Probability: frequentist and Bayesian definitions

What is a probability? You might be surprised to learn that statisticians don’t really know. In the frequentist philosophy of statistics², a probability is the fraction of times a certain outcome would be expected in the long-run. So, for example, $p(heads)=0.5$ when you flip a coin because we would expect to get 50 heads if we flipped a coin 100 times. The frequentist definition depends on the ability to do something multiple times, and is in no way speculative. Thus a frequentist can’t give a probability of life on Mars, because we haven’t been to any other planets and so a frequentist has no idea what the overall ‘rate’ of life being found on planets is. Frequentists are so anally-retentive about this that one professor, finding themself imprisonned during World War II, decided to toss a coin ten thousand times to make sure there are equal probabilities of getting a head or a tail³.

Under the Bayesian philosophy, a probability is a belief about the relative likelihood of something happening, which can be modified in the face of data. If that seems a bit hand-wave-y, that’s because it is: it’s a belief, and so it has no connection with any formal quantitative number. It’s possible, in the Bayesian framework, to make relative comparisons between the probabilities of things (“this is twice as likely as that”), but ultimately there’s no connection to a hard, externally-verifiable reality. While this might sound a bit terrifying, it shouldn’t be surprising: the frequentist definition is rooted in empiricism, and so that base is already covered. This may all seem so strange that there seems no advantage to a Bayesian perspective; there is, and it all relates to the concept of relative probability, but we’ve not got time to go into it yet. Just be aware that there are, fundamentally, only⁴ two kinds of probability: one rooted in the relative frequency of events, and the other in belief⁵.

1.2 Statistical distributions

A statistical distribution is easier to intuitively understand from a frequentist perspective. A statistical distribution describes the relative frequencies with which samples of given values will be drawn. Imagine we were to go out and record (sample) the heights (values) of 100 undergraduates and plot them in a histogram. This histogram would have a particular shape, and we could describe that shape using a mathematical equation. If we divided the vertical axis by the total number of samples collected (leaving the shape of the histogram unchanged), we would have a plot of the distribution of the relative frequencies of each height. A frequentist would say that, the more heights you sampled, the closer that distribution would get to the true statistical distribution of heights.

Statisticians often say that samples are drawn from a distribution. This means that they find it convenient to statistically model a given variable (e.g., height) as if all its samples follow the distribution they specify (Figure 1.1). How statisticians go about drawing new samples from a distribution is a complicated business. Frequentists often use associations between physical objects and statistical mechanics: in other words, they do carry out some physical process that obeys the same mathematical principles as their statistics and simultaneously introduces ‘randomness’ to draw from their distribution of choice. A good example of this is tossing a coin, another is the throwing of beads into cups⁶. We are lucky that modern computers have ‘randomness’ in the form of electrical fluctuations that we can use to draw numbers from a distribution⁷.

Some distributions are continuous, such as the height example above, such that it only makes sense to think of ranges within them. For example, we might ask informally what fraction of the population we might expect to have a height of 190 cm, but really we’d know that’s a shorthand for asking what fraction of the population we’d expect to have a height greater than or equal to 1.85 meters and less than 195 cm. If we wanted to know the fraction of our sample that is 190cm tall, we would count the number of people between 185 and 195 cm tall and divide that by the total number of people we sampled. We could also read that off our histogram: (1) find the box(es) that cover people between 185 and 195 cm tall, (2) multiply their heights (the frequency) by their widths to get the total number of people in that range, (3) add up the boxes, and then finally (4) divide by the total number of people surveyed. The normal distribution is continuous. Some distributions are discrete and don’t display this property; there are many examples, but the simplest is maybe the Poisson distribution that is only capable of spewing out whole numbers.

The problem with statistical distributions is they have a tendency to look terrifyingly complex. I think it’s important to remember that each wasn’t ‘made up’ to look a particular way, but rather was mathematically derived to describe samples that were taken under certain conditions. The most important example of this is the normal distribution, which was derived by the ‘Prince of Mathematicians’ Gauss⁹. Gauss showed that when you sample something that is affected by lots of independent errors, the distribution of those samples ($x$) will be centred around the true mean ($\mu$¹⁰) with a spread defined by the standard deviation ($\sigma$¹¹) as in this equation:

\[ f(x) = \frac{1}{\sqrt{2 \pi} \sigma} e ^\frac{-(x-\mu)^2}{2 \sigma^2} \tag{1.1}\]

This equation defines the normal distribution, and the theorem that distributions will converge on it like this is central to modern statistics and is called the Central Limit Theorem. The equation is complicated, and you don’t need to learn it, but it’s important that you remember that this is why so many histograms tend to look bell-shaped (e.g., Figure 1.1). Things like people’s heights are affected by thousands of different factors (e.g., hundreds of genes, dietary choices, lifestyle, and more), and so all those differences add up, through the Central Limit Theorem, to make people’s heights normally distributed. Statisticians get so excited about the normal distribution precisely because the Central Limit Theorem points to it: otherwise, it would just be a bizarrely complicated way of drawing a reasonably pretty curve.

In Equation 1.1 above, I referred to a normal distribution’s parameters: its mean ($\mu$) and standard deviation ($\sigma$). The mean is the centre of the distribution, and most of the samples from that distribution will look more-or-less like that value (it’s the distribution’s ‘central tendency’ that, in the ‘limit’ where you sample from it many times, it will centre on). The standard deviation is the spread of the distribution, and observations from a normal distribution with a greater value will exhibit more variation. You’ve probably encountered these terms before, but it’s critical that you understand there’s nothing magical about them. They just happen to be the parameters of the normal distribution; if we were modelling data from another distribution we might use totally different parameters. Indeed, some distributions don’t even have defined means or standard deviations!

Because frequentists believe that true distributions and probabilities come from sampling, they don’t really need maths to find the normal distribution (or any distribution, for that matter). Frequentists believe that, if they keep sampling, the distribution of values they get will converge on the true distribution. Take a look at Figure 1.1, where I show you how much nicer (smoother) and more like the true, underlying distribution, the distribution of samples come when I sample more. If I sampled even more (the grey boxes in Figure 1.2), I could make the width of the histogram bars narrower: I’ve got more data, so I can split it more finely, and it more accurately shows the shape of the distribution. You would probably agree with me that I could sample even more times, and make the bins even narrower, with the only caveat being I’d have to change the scale on my vertical axis (look at the heights on Figure 1.1; I have to keep changing the limits as I sample more). To derive distributoins, frequentists write a mathematical equation where they sample an infinite number of times. This allows them to make the bars infinitely narrow (because they’re doing loads of sampling), and allows them to draw a curve along the mid-points of those bars (they’re infinitely narrow, so really they form a curve anyway).

So how do we figure out the probability of a particular observation? Remember how height isn’t really a single number? How when we thought about people 190 cm tall, we really were asking how many were taller than 185 cm and shorter than 195 cm? It’s exacxtly the same with our statistical distributions, becuause to a frequentist the underlying logic and maths is all the same. If we want to know the probability of someone being 190 cm tall, we can read that off our histogram: (1) find the box(es) that cover people between 1.85 and 1.95 meters tall, (2) multiply their heights (the frequency) by their widths to get the total number of people in that range, (3) add up the boxes, and then finally (4) divide by the total number of people surveyed. The only difference now is that the boxes are infinitely narrow, and so it’s easier to draw a curve rather than drawing loads of boxes (Figure 1.2). When plotting our curve, just as we had to rescale our histograms, now we have to rescale for infinite sampling, and so divide our frequencies so that the total area under the curve is exactly equal to 1. We need that total area under the curve to be 1 because it’s a probability distribution: we need a 100% probability that something will happen, and everything that can happen needs to be on the curve. We call this plot a probability density function because we use them to get probabilities, just as in histograms we are interested in the area under the curve (that’s the ‘density’ part), and it’s defined by an equation we’ve derived assuming infinite sampling so we can call it a ‘curve’. The probability density function for the normal distribution is given to you in @equ-normal. Remember: a probability is just the area under the curve.

Figure 1.2: **Frequentists derive true statistical distributions by sampling towards infinity.** A frequentist can derive a statistical distribution by deriving, mathematically, what happens when you sample an infinite numbers of times. Our histogram’s bins become infinitely narrow, and we can draw a curve across the top of them (the red line at the right). We can calculate the probability of any set of observations (*e.g.*, a height of 190 cm) by measuring the area under the curve between two points that define your question (*e.g.*, the area between 185 and 195 cm).

1.3 Estimating parameters: maximum likelihood

So far we’ve talked about drawing samples from a distribution with known parameters. You may be wondering how, given a set of samples, we can determine the distribution (and parameters) most likely to have generated those samples. We’ve already talked about estimating the probability of observing a range of possible values (e.g., a height between 185 and 195 cm) and, as outlined above, we can’t get a probability of an exact height because it’s a continuous variable. We can, however, read off probability density of an exact number (e.g., 190 cm): that’s the likelihood. For example, plug in an observation ($x$) and parameters ($\mu$ and $\sigma$) to Equation 1.1, and voilà, you’ve got the likelihood that your observation was drawn from that normal distribution. Thus to find out what distribution/parameters a set of samples were drawn from, we find the distribution/parameters that maximise the product of the likelihoods of all the datapoints. We take the product (‘product’ means ‘multiply’) because it’s the probability we got each value together, not independently, and the maximum because bigger means more likely/probably. We call this the maximum likelihood solution, or ML estimate of a distribution. Notice that it’s an estimate, taken from the sample, of the parameter of the true distribution from which those samples were drawn. The true value of the underlying, true distribution is called the population estimate.

Estimating all these parameters has a cost. Intuitively, you would probably agree that the more things we estimate, the more underlying data we would need. In statistics, the concept of degrees of freedom helps us keep track of this. Imagine a friend sampled ten heights from a normal distribution, told you nine of them, and then offered you a million pounds if you exactly guessed the tenth. You’d say they were being unfair, because there’s no way you could get it exactly right—you could get close, but not exactly there. Now imagine that they gave you the mean of all ten numbers as well: you would laugh, because using simple algebra¹² you could exactly predict the missing value using the formula for a mean¹³. The total degrees of freedom of a random sample is equal to the number of data points in it; once you estimate a parameter from that sample, you ‘soak up’ (use) a degree of freedom and so you can predict one of those data points from all the others and that estimate. In the previous example, calculating the mean of the same meant there were 9 remaining degrees of freedom ($10 - 1$; 10 data points minus one estimated parameter). Thus you always need at least one data point for each parameter you estimate, and Very Bad Things will happen if you don’t have that. My general piece of advice is to make sure you have at least ten data points per estimated parameter, but that’s based less on statistical rigour and more on experience.

The problem with a single ML solution/estimate is it doesn’t give you any idea of the uncertainty associated with that estimate. One way to get at this is with confidence intervals, often abbreviated to CIs. Confidence intervals allow us to bracket the uncertainty around a parameter estimate that we’ve found through an ML search. To find them, you have to make a decision: how worried are you about making an error? If you want to be really, really certain that you’re going to construct some confidence intervals that will contain the true parameter value, then you want to have very wide CIs. The problem is, you won’t then have a lot of power to know with any certainty what the true value is: the interval will be so wide that it will encompass a lot of possible values. Make sure you take a look at Figure 1.3, which gives you an example of what all of this looks like; we’ll go through this in more detail in the next chapter.

In general, the more data you have, the greater your power to detect something that’s going on in your data. The more complex your model—the more parameters that you need to build your statistical distribution—the more degrees of freedom you will soak up, and so the more data you will need. A statistical model that misses a real effect (e.g., mood is affected by time of day, but you can’t find evidence for it in your data) has made a Type II error. Conversely, if you find evidence for an effect that isn’t really there (e.g., you suggest there is a link between the MMR vaccine and autism, but really there isn’t) you have committed the opposite kind of mistake: a Type I error. The way to remember the difference is the story of ‘The Boy Who Cried Wolf’: first he cries there’s a wolf when there isn’t one (Type I error), then he cries wolf when there is and no one believes him (they commit a Type II error).

Figure 1.3: **A graphical overview of Maximum Likelihood estimation.** As described in the text, it’s possible to find the best estimate of a given parameter by calculating the likelihood of the observed data if it were drawn from a distribution parameterised by different estimates. Doing so creates a *likelihood curve*, whose *maximum* represents the best-estimate of that particular parameter (the Maximum Likelihood Estimate, or MLE). It is also possible to construct *confidence intervals* by going a certain distance away from that maximum. Note that, because likelihoods are almost always incredibly small, we work with them on the log scale. Note that it is not coincidence that this likelihood curve looks approximately normal; if you should ever look at an empirical ML curve and it does not, be careful, as it likely indicates something has gone Badly Wrong.

1.4 The zoo of statistical distributions

When I was a graduate student, I used to almost ‘collect’ statistical distributions, and felt that the more I knew, the ‘better’ I was at statistics. Sadly, I was wrong¹⁴: it’s better to understand a few distributions, and where they come from, than to collect them like imaginary pocket monsters. With that in mind, I’m going to give you a brief tour of some useful statistical distributions, where they come from, and what they can be used for. I’m going to assume that you’re comfortable with the concept of a normal distribution, either informally (from the text above) or formally (from the extension section below where we derive it from first principles). I’m also not going to cover the $t$ distribution because it’s the basis of the $t$-test, which is what we will be covering extensively in the next chapter. Note also that I’m not giving you nice pictures of what these distributions look like: doing so is one of your exercises below…

1.4.1 Chi-squared ($\chi^2$)

A $\chi^2$ distribution is what you expect to get if you take the variances of normal distributions and sum them up. Let’s state that semi-formally¹⁵:

\[ Q = \sum \sum (x_k - \bar{x_k})^2 \tag{1.2}\]

Where Q is a $\chi^2$-distributed variable, $x_k$ is an observation of the $k$th random normally-distributed variable we’re working with, and $\bar{x_k}$ is the mean of the $k$th normal distribution. Note that there are two summation signs here because we’re doing two sets of summing¹⁶: once within each random variable, summing the squared deviation from the mean (the variance) for each draw within that sample, and then again to sum across all the variables’ squared differences. To calculate the probability density function for this variable, you would carry out that operation across the probability density function of a set of normal distributions.

It’s a useful distribution to know about because it comes up so much; it’s essentially a measure of goodness of fit. Now you know how it’s derived from the variances, perhaps that makes some sense, or at least perhaps it will after the next chapter where we will be calculating variances over and over again until you’re sick of them.

1.4.2 F

The F-distribution is so-named for Ronald Fisher, who worked extensively on model-testing, and it’s incredibly important in model testing and we will spend a lot of time working with it next session. Sadly, Fisher was a horrible person, and notably racist even in a time when racism was rampant and widely accepted, so I don’t want to dwell on the distribution’s origins anymore. Luckily, the distribution itself is quite simple¹⁷, albeit looking a bit weird because it is two fractions ($\frac{Q_1}{d_2}$ and $\frac{Q_2}{d_2}$) that are themselves in a larger fraction:

\[ F = \frac{\frac{Q_1}{d_1}}{\frac{Q_2}{d_2}} = \frac{Q_1d_2}{Q_2d_1} \tag{1.3}\]

Where $Q_1$ and $Q_2$ are two $\chi^2$-distributed variables (see above) and $d_1$ and $d_2$ are the degrees of freedom of those variables. As we’ll see next time, this is a very useful distribution because it tells you whether there’s anything unusual about the differences in variation of two distributions. It’s unusual to see it written out in terms of two sets of fractions, but I encourage you to think of it like that as it helps simplify the concepts we will encounter next time.

1.4.3 The Poisson distribution

The Poisson distribution is the first ‘oh wow’ distribution for many, because it solves a problem you have been repeatedly told, as an undergraduate, wasn’t a problem but in reality actually is a Very Big Deal: nearly-but-not-quite-continuous variables. Not everything in life can take non-integer values: lightning can strike once, twice, thirty times, but not half a time. The Poisson distribution can deal with this, and is predicts only whole numbers:

\[ p(x=k) = \frac{\lambda^ke^{-\lambda}}{k!} \tag{1.4}\]

What this equation means is that the probability that you observe $k$ things (lightning strikes, horse kicks, fish swimming past you) in your random variable $x$ is a function of the natural number $e$ and $\lambda$. $\lambda$¹⁸ is the only parameter in the Poisson distribution, and is equal to both the mean and the variance of the distribution. We need all the fancy $p(x=k)$ on the left-hand side because we must calculate discrete probabilities for each value of $k$, and there can be no $k$ of 1.5 (only 1 or 2), and thus the Poisson distribution is not really technically a continuous distribution. This all makes the maths of deriving the distribution a bit involved, so we won’t dwell on it.

1.4.4 The Binomial distribution

Finally, the workhorse of the bioinformatic world: the Binomial distribution. It’s the probability of seeing an event happen a certain number of times given how many times you looked and the underlying probability of it happening. For example, it can give you the probability of getting ten heads when you toss twenty coins. It can give you the probability of seeing 5 diseased individuals in a population of 100 given a background prevlance of 1%. It’s your best friend, really it is, and it’s so commonly used it’s got its own special notation:

\[ p(x=k) = \binom{n}{k} p^k (1-p)^{n-k} \tag{1.5}\]

… where …

\[ \binom{n}{k} = \frac{n!}{k!(n-k)!} \tag{1.6}\]

…I know, right? This is definitely the worst one to look at the maths of. Bear with me. It’s giving us the probability of detecting $k$ ‘successes’, where a success is ‘the event we’re interested in happened’ and is quite often, in my experience, something that doesn’t feel like a success such as ‘someone has a disease’. We are doing the thing (testing a patient, tossing a coin) $n$ times, and our probability of success is $p$. Some people introduce the term $q$ in the above equations to represent the probability of a failure, but I think not writing it out makes the logic of Equation 1.5 easier to track through: product of the probability you succeed, the probability you fail (which is just when you don’t succeed), and a correction for how many ways there are of doing the two. Let’s do it with a worked example: tossing four coins and getting two heads. It’s the probability of all the success you would need ($p^k$; two heads are $0.5 \times 0.5 = 0.5^2 = 0.25$) multiplied by the failures you would need ($(1-p)^{n-k}$; two tails are $(1-0.5) \times (1-0.5) = (1-0.5)^2 = 0.5^2 = 0.25$). We then correct for all all the different ways of throwing two heads in a row (where H is ‘heads’ and T is ‘tails’: HHTT, HTHT, HTTH, TTHH, THTH, THHT—six ways total): $\frac{n!}{k!(n-k)!} = \frac{4!}{2!(4-2)!} = \frac{4!}{2!2!} = \frac{4 \times 3 \times 2 \times 1}{(2 \times 1) \times (2 \times 1)} = \frac{24}{4} = 6$. In total: $0.25 \times 0.25 \times 6 = 0.375$.

Look, it’s a lot of maths. But trust me, it’s a lot easier than drawing out those blasted branching diagrams you may have been made to do previously at A-level… Oh yes and, of course, will be made to do in the exercises. Sorry about that.

1.5 Extension: Deriving the normal distribution

There is no need to read this material, particularly on a first read or if you’re not very comfortable with mathematics. If you are, however, you may find this interesting. There are plenty of derivations of the normal distribution online, and if you are interested (and have the mathematical background) I would encourage you to look a couple up. It is possible, however, with the equivalent of an A-level in maths, to understand the important principles of the derivation of the normal distribution and I think it is a useful exercise to go through because it somewhat demystifies what the normal distribution represents. It didn’t fall from the sky and it wasn’t made up for fun, it follows from relatively simple assumptions. With that in mind, what follows is an attempt to capture the spirit of a full derivation.

To say a variable is normally distributed is to say that the frequency of observations drawn from that distribution is lesser the further away from the mean of those observations. Let’s write that in maths:

\[ \frac{df}{dx} = -v(x - \mu) f(x) \tag{1.7}\]

Where $f(x)$ is the normal distribution we’re trying to derive, $\mu$ is the mean, and $v$¹⁹ is a measure of how fast the frequency falls off away from that mean. The left-hand side of the equation is a derivative: we’re talking about the rate of change of the distribution as we move along the values of the observation ($x$), so we’re starting from this change (our starting observation that motivated the distribution in the first place) to calculate the actual distribution itself ($f(x)$). There’s $x-\mu$ because we are defining the shape of the distribution by how far away we are ($x$) from the mean ($\mu$), and we multiply that by $-v$ because as we ‘go away from’ the mean (become more positive than the mean) we want the frequency to go down (and go up as we approach it—become less negative). Finally, we multiply everything by $f(x)$ because we want everything to be proportional to the frequency at that point along the distribution: if we have a very low probability, we want to change very slowly because we never want to hit zero and have impossible values, and we want to change fast closer to the centre of the distribution. Equally, if you’ve been paying attention earlier, you remember that probability distributions need to have an area equal to 1, and so having this in here is convenient in that sense too (but more on that later).

OK, I think the above was the most challenging part so if you’re still with me well done. Now let’s re-arrange a bit to get all our $x$s on one side:

\[ \frac{df}{f} = -v(x-\mu) dx \tag{1.8}\]

That wasn’t so bad! We just multiplied by $x$²⁰! Now let’s get rid of those nasty derivatives by integrating both sides:

\[ \int \frac{df}{f} = -v \int (x-\mu) dx \tag{1.9}\]

The integral on the left-hand side is something you can look up in a table (or derive from first principles if you’re that good at maths but that bad at remembering things), and the right-hand side is not so bad if you’ve done calculus but I do appreciate it’s not trivial. Once we’ve done that we have:

\[ \ln f = \frac{-v(x-\mu)^2}{2} + c \tag{1.10}\]

…where $c$ is an integration constant as per usual. Let’s exponentiate to finish this off and finally write this in terms of $f(x)$:

\[ f(x) = c e^{\frac{v}{2}(x-\mu)^2} \tag{1.11}\]

Great! Now this nearly looks like Equation 1.1, but it isn’t quite there, is it? In particular, we have that nasty integration constant $c$ lying around, and we also have $v$. Let’s tackle $c$ first, because it’s mathematically the hardest part to deal with (and is where we’ll stop with the maths) but is conceptually the simplest.

Remember that all probability distributions have to have an area under them equal to 1? This is because something has to happen, and so the total realm of possible probabilities has to sum to one. You would probably agree with me that Equation 1.11 cannot have a negative value (there is an $e$ in there) unless $c$ is negative, thus $c$ is essentially a normalising constant that determines what the area under that curve will be. Thus the tricky part of the maths of the derivation of the normal distribution is finding a general solution for $c$ such that the integral of Equation 1.11 (i.e., the area under the curve) always sums to 1 for any given value of $\mu$ or $v$. That’s all there is to it; it happens that the solution is to set $c$ to be $\frac{1}{\sqrt{2\pi\sigma}}$ (I’ll get to what $\sigma$ is…) but there’s no great magic to it, just a lot of fiddling. So now you know where that normalising constant comes from and what it achieves.

So now to $v$. Remember that $v$ is a way of measuring how quickly we swoop up/down when we come closer to or further away from the mean ($\mu$). I suppose, therefore, it a measure of how variable the distribution is: larger values of $v$ will mean that the distribution moves more ‘quickly’ around the mean. Hmmm, so I guess that’s a bit silly to call it a measure of how variable the distribution is, then, because bigger values of $v$ will tend to be in distributions that are less variable (faster sweeps up/down). Why don’t we fix that—it’s just a variable, after all—and replace it with its inverse ($\frac{1}{v}$)? That would give us:

\[ f(x) = c e^{\frac{1}{2v}(x-\mu)^2} \tag{1.12}\]

Great! Oh, hang on, ‘variable’ isn’t a good word to use because we’re studying random variables and that’ll be really confusing. Let’s call it variance instead. Oh, hang on, we can’t use $v$ because we use Greek letters for parameters (like the mean—$\mu$) in statistical distributions because otherwise the maths looks much too approachable and we can’t scare graduate students. Also, $v$ looks like the Greek letter $\nu$, which is actually said as an ‘n’ sound. Oh, stuff it, let’s go for something silly—let’s call it $\sigma$. Oh wait, hang on, we can never have that variable be negative²¹ because otherwise we’ll come ‘down’ towards the mean and then go back ‘up’, and have a distribution that is weirdly shaped. I guess an easy way to ensure that never happens is just to square it, right, because squared values can never be negative²², right? OK, let’s call it $\sigma^2$ then:

\[ f(x) = c e^{\frac{1}{2\sigma^2}(x-\mu)^2} \tag{1.13}\]

Oh look! We’re done! Hooray!

1.6 Exercises

Today’s exercises are going to seem a little strange, because we don’t know enough practical statistics to start working with real datasets yet. Focus on getting comfortable with the concepts, and use these questions as ways to test your understanding. These are the most difficult questions you will face in the course, so don’t worry if they seem tough.

Fundamentals of probability definitions. If you’re writing answers to these, do not write pages and pages for these questions: you should be able to answer them in a few sentences.
1. Your friend Sarah warns you that there’s a 10% chance (in a frequentist sense) that Mount St Helens will erupt tomorrow. Your other friend, Michael, scoffs, and says that her probability makes no sense because the volcano will either erupt tomorrow or it won’t, and because there is only one tomorrow there can be no frequentist probability of that. What do you think?
2. Your friend Michael approaches you with a problem. He has the average blood pressures of two professors, and wants to statistically test whether the professors’ pressures are drawn from a Blaargian or a Floogian distribution. The probabilities each measurement are from a Blaargian distribution are $0.1$ and $0.9$, and the probabilities the pressures are from a Floogian distribution are $0.4$ and $0.6$. Which distribution do you think the data are most likely drawn from? Show your working (hint: if the working isn’t trivial, you’ve gone wrong).
3. Michael is beginning to disturb you a little, and wants to keep talking about these blood pressures. He tells you that the Blaargian distribution is defined by one parameter ($\alpha$), while the Floogian distribution is defined by two parameters ($\alpha$ and $\beta$). Does this affect your answer to question 1b above? Why?
Going to the statistical zoo. The aim of this exercise is not to produce hundreds of plots that you laborious copy-paste into Word, but rather to become familiar with both some of R’s basic plotting functions and also the shapes of statistical distributions. Use the code below to generate draws from statistical distributions (which you can lookup in R help with something like ?rnorm), and then sketch the impact of changing different parameters either on a piece of paper or, if it is easier for you, using line drawings in Word or something like that. Make sure you look carefully at the horizontal (the ‘x’ axis) in these plots, because otherwise you will miss differences between plots.
1. The normal distribution has two parameters: $\mu$ and $\sigma$. Use the code below to investigate the impact of changing $\mu$ on the normal distribution, and sketch out your findings.
  # Draw 10,000 times from a normal distribution with # a mean (mu) of 0 and sd (sigma) of 1 # and then make a histogram of those samples hist(rnorm(n=10000, mean=0, sd=1))
  # Plot out again, changing the mean hist(rnorm(n=10000, mean=5, sd=1))
  # Change the plot so that we can have two plots on the screen # at the same time, then plot the histograms side-by-side par(mfrow=c(1,2)) hist(rnorm(n=10000, mean=0, sd=1), main="mean=0") hist(rnorm(n=10000, mean=5, sd=1), main="mean=5")
  # Reset the plotting window (clearing the two plots) # so that you can plot one thing at a time again par(mfrow=c(1,1))
2. Use the code above to sketch the impact of changing $\sigma$ (note that R uses the standard deviation, $\sigma$, to parameterise the normal and not the variance, $\sigma^2$ but, as one is just the square of the other, it needn’t bother you).
3. The Poisson distirbution has just one parameter, $\lambda$. It can be plotted out like this:
  hist(rpois(n=10000, lambda=1))
  …sketch out how it changes across $\lambda$. Pay attention to the asymmetry of its distribution.
4. The $\chi^2$ distribution is parameterised by the number of distributions it’s summing over (there are caveats on this; see the next few chapters). Use the R function rchisq, and its two arguments n and df, to explore how the distribution changes across values of df such as 1, 5, 10, 50, and 100. Note, again, its asymmetry. Are you beginning to see a pattern in how R defines its functions that draw Random numbers from a distribution?
5. Drawing random numbers from the Binomial distribution often confuses students because it requires three arguments: the number of times to run the simulation (draw numbers) (n), the numbers of times to attempt a success within each draw/run (size), and the probability of success (prob). Make a sketch of how different values of size affect how ‘discrete’ and ‘gappy’ the distribution looks, but first follow the outline of how to explore the function given below.
  # Let's toss a coin in the air ten times rbinom(n=10, size=1, prob=0.5) #> [1] 1 0 0 1 0 1 1 1 1 1 # Let's toss ten coins in the air and count the number of heads rbinom(n=1, size=10, prob=0.5) #> [1] 7 # Let's toss ten coins in the air and count the number of heads ten times rbinom(n=10, size=10, prob=.5) #> [1] 6 6 4 3 7 4 3 5 4 6 # ... make sense? Great. If not... Play around a bit, then ask me.
6. The F distribution has too many parameters to be easy to sketch, so take five on that one. This question is left deliberately blank.
The Binomial distribution. You cannot become familiar enough with the Binomial distribution; it really does make your life much easier. While it is always possible to find precise, analytical solutions to Binomial problems, you may find that learning a thing or two about how to simulate draws from it comes in handy from time to time.
1. Sarah hands you a bag with five marbles in it, two of which are red and the other three black. She says you can pull a marble out and put it back in the bag (i.e., in statistical parlance you are ‘drawing with replacement’) three times. She asks to calculate, without using the Binomial distribution but by drawing a probability tree diagram, the probability you will draw three red balls.
2. She asks you to now estimate the probability you will draw two red balls, again only using a tree diagram. You are so glad that she helped you with that whole blood pressure situation with Michael that, honestly, you just go along with it.
3. Sarah, who makes it clear to you that she finds tree diagrams boring as well, asks you to use the Binomial probability equation your biostatistics lecturer gave you to estimate the answers to (a) and (b) for her. She says her lecturer never gave it to her and she feels hard-done-by, and so would like you to do her a favour.
4. After calculating the probabilities by hand, Sarah is not so sure she was that hard-done-by. Sarah shares a trick her lecturer taught her with you, which is to use the Binomial distribution density function (dbinom) to estimate the probabilities from parts (a) and (b). She asks if you can do it and, not wanting to look foolish, you glance back at your workings from question 2e and give it a go.
5. Sarah’s so impressed with your ability to think on your feet, she asks if you can work together to estimate the probability of drawing (with replacement) between 30–50 red balls (inclusive) from a bag with 200 balls in total, of which 30 are red and the remaining 170 are black. She knows lots of ways to do it, but feels there’s got to be an easy ‘one-liner’ in R that could do it for her. Hint: Take a look in the help entry for dbinom (?dbinom).
6. Sarah notices a question posted by Michael (who you have both now decided to aviod)written in the margin of her class workbook, wondering whether it matters if half of the black balls in the bag are replaced with blue balls. Would this affect any of your answers? Why (not)?
7. Sarah notices a final entry in her workbook from Michael, which appears to be somewhat hastily written and followed by frantic scratchings. It seems to follow a question about whether it would be possible to repeat these exercises ‘without replacement’: in cases where no ball taken from the bag was returned afterwards. Draw out the probability diagram for part (a) and, on the basis of that, consider whether it was wise of Michael to attempt part (e) using the Binomial distribution.

Easier said than done, unless of course you’ve never been told anything about statistics…↩︎
Note that no one statistician is ‘a frequentist’ or ‘a Bayesian’, regardless of what you might read online. To describe yourself as ‘a frequentist’ is a bit like saying you’re ‘a portrait artist’. That’s great and everything, but if portraits are literally all you can paint then you’re probably not much of an artist. It is, however, a useful shorthand when writing handouts, so I’m going to use it occasionally. Just don’t introduce yourself to a real statistician as ‘a Bayesian’ to try and sound clever, because I can assure you that you won’t.↩︎
And what, exactly, did you do during the pandemic, eh? Kerrich JE (1946; An experimental introduction to the theory of probability, E. Munksgaard). Believe it or not, people are still arguing about this! See Gelman & Nolan (2002; The American Statistician 56(4):308).↩︎
…there are a arguably a few more, but let’s not worry about that right now…↩︎
Bayesian philosophy is named for Reverend Bayes—it shouldn’t be surprising that a priest invented a statistical method based around belief…↩︎
Several mathematicians have ‘proven’ theorems by throwing beads into cups; look up the “Galton board” (sometimes called a ‘quincunx’). If there ever were a reason to be grateful for calculus, this is that.↩︎
The commonest ways of drawing random numbers (including the ones used by R) are done using a different approach, and you ask me if you are interested.↩︎
It was a very busy weekend, let me tell you↩︎
And so the normal distribution is often called the Gaussian distribution; I will use ‘normal’ because I think it sounds like scary..↩︎
‘mu’—the Greek letter for ‘m’ (mean)↩︎
‘sigma’—the Greek letter for ‘s’↩︎
Perhaps not simple, but for a million pounds?↩︎
…$\frac{\sum{x}}{n}$…↩︎
…as I found out in exams↩︎
That phrase ‘semi-formally’ is doing quite a lot of work here. This is less a definition and more the start of a derivation, to be completely honest, but if you’ve noticed that then you don’t really need to be reading this part of the book anyway…↩︎
…I’m not putting in the formally required subscripts on the summation signs for simplicity↩︎
Ahem. Well, it is when you do a simplification along the lines of the trick I used above… This is the last trick in this section, I promise!↩︎
‘lambda’—the Greek letter for ‘l’↩︎
I’m using non-standard terminology for this constant but you will see why at the end.↩︎
Alright, I know calculus doesn’t quite work like that, but it really does look that way…↩︎
This comes from the maths that I ignore about deriving $c$, but the argument I’m giving here is valid↩︎
Shut up at the back, I know they sometimes can, but we’re dealing with real random numbers here, OK? Give me a break…↩︎

# Philosophy of statistics In this first chapter, I'm going to try and give you an overview of the core concepts that define statistics. The only way this can possibly work is if you bear with me and promise to try and forget everything you think you have been told about statistics^[Easier said than done, unless of course you've never been told anything about statistics...]. It's easy to get the wrong end of the stick with these concepts, particularly if you first encountered them some time ago, and often the easiest way to learn is to start again. If you're coming from a discipline other than statistics, please also try and bear in mind that statistics is likely a larger discipline than your own subfield. This means that, in order to develop a serious statistical intuition and understanding, you've got to learn some fundamental vocabulary and concepts. Don't approach this class with the mindset of trying to figure out how it will help you analyse the dataset you have in front of you right now: treat it as a grounding in the fundamental concepts that will allow you to learn anything you want about any dataset. If you're impatient and just want answers, then fair enough, but please just bear with me... ## Probability: frequentist and Bayesian definitions What is a probability? You might be surprised to learn that statisticians don't really know. In the *frequentist* philosophy of statistics^[Note that no one statistician is 'a frequentist' or 'a Bayesian', regardless of what you might read online. To describe yourself as 'a frequentist' is a bit like saying you're 'a portrait artist'. That's great and everything, but if portraits are literally all you can paint then you're probably not much of an artist. It is, however, a useful shorthand when writing handouts, so I'm going to use it occasionally. Just don't introduce yourself to a real statistician as 'a Bayesian' to try and sound clever, because I can assure you that you won't.], a probability is the *fraction of times a certain outcome would be expected in the long-run*. So, for example, $p(heads)=0.5$ when you flip a coin because we would expect to get 50 heads if we flipped a coin 100 times. The frequentist definition depends on the ability to do something multiple times, and is in no way speculative. Thus a frequentist can't give a probability of life on Mars, because we haven't been to any other planets and so a frequentist has no idea what the overall 'rate' of life being found on planets is. Frequentists are so anally-retentive about this that one professor, finding themself imprisonned during World War II, decided to toss a coin ten thousand times to make sure there are equal probabilities of getting a head or a tail^[And what, exactly, did you do during the pandemic, eh? Kerrich JE (1946; *An experimental introduction to the theory of probability*, E. Munksgaard). Believe it or not, people are still arguing about this! See [Gelman & Nolan (2002; The American Statistician 56(4):308)](https://www.tandfonline.com/doi/abs/10.1198/000313002605).]. Under the *Bayesian* philosophy, a probability is a *belief* about the relative likelihood of something happening, which can be modified in the face of data. If that seems a bit hand-wave-y, that's because it is: it's a belief, and so it has no connection with any formal quantitative number. It's possible, in the Bayesian framework, to make relative comparisons between the probabilities of things ("this is twice as likely as that"), but ultimately there's no connection to a hard, externally-verifiable reality. While this might sound a bit terrifying, it shouldn't be surprising: the frequentist definition is rooted in empiricism, and so that base is already covered. This may all seem so strange that there seems no advantage to a Bayesian perspective; there is, and it all relates to the concept of *relative* probability, but we've not got time to go into it yet. Just be aware that there are, fundamentally, only^[...there are a arguably a few more, but let's not worry about that right now...] two kinds of probability: one rooted in the relative *frequency* of events, and the other in belief^[Bayesian philosophy is named for Reverend Bayes—it shouldn't be surprising that a priest invented a statistical method based around belief...]. ## Statistical distributions A statistical distribution is easier to intuitively understand from a frequentist perspective. *A statistical distribution describes the relative frequencies with which samples of given values will be drawn.* Imagine we were to go out and record (*sample*) the heights (*values*) of 100 undergraduates and plot them in a histogram. This histogram would have a particular shape, and we could describe that shape using a mathematical equation. If we divided the vertical axis by the total number of samples collected (leaving the shape of the histogram unchanged), we would have a plot of the distribution of the relative frequencies of each height. A frequentist would say that, the more heights you sampled, the closer that distribution would get to the true statistical distribution of heights. Statisticians often say that samples are *drawn* from a distribution. This means that they find it convenient to statistically model a given variable (*e.g.*, height) as if all its samples follow the distribution they specify (@fig-sampling). How statisticians go about drawing new samples from a distribution is a complicated business. Frequentists often use associations between physical objects and statistical mechanics: in other words, they do carry out some physical process that obeys the same mathematical principles as their statistics and simultaneously introduces 'randomness' to draw from their distribution of choice. A good example of this is tossing a coin, another is the throwing of beads into cups^[Several mathematicians have 'proven' theorems by throwing beads into cups; look up the "Galton board" (sometimes called a 'quincunx'). If there ever were a reason to be grateful for calculus, this is that.]. We are lucky that modern computers have 'randomness' in the form of electrical fluctuations that we can use to draw numbers from a distribution^[The commonest ways of drawing random numbers (including the ones used by `R`) are done using a different approach, and you ask me if you are interested.]. ::: {#fig-sampling layout-ncol=2} ![Regression](Figs/01-hist-small){#fig-sampling-small} ![ANOVA](Figs/01-hist-bigger){#fig-sampling-small} **Frequentists believe that when you sample more, your data's distribution approaches the true distribution.** In the histogram on the left, I've measured 50 people's heights and plotted a histogram of the data. In the histogram on the right, I've sampled one million people's heights^[It was a very busy weekend, let me tell you] and plotted their histogram. A frequentist believes that the more samples you take, the closer your distribution gets to the distribution of the total population you're sampling from. That's why the smaller sample on the left looks much less believable than the one on the right (it has that weird peak, for example): the more we sample, the more like reality everything looks. ::: Some distributions are *continuous*, such as the height example above, such that it only makes sense to think of ranges within them. For example, we might ask informally what fraction of the population we might expect to have a height of 190 cm, but really we'd know that's a shorthand for asking what fraction of the population we'd expect to have a height greater than or equal to 1.85 meters and less than 195 cm. If we wanted to know the fraction of our sample that is 190cm tall, we would count the number of people between 185 and 195 cm tall and divide that by the total number of people we sampled. We could also read that off our histogram: (1) find the box(es) that cover people between 185 and 195 cm tall, (2) multiply their heights (the frequency) by their widths to get the total number of people in that range, (3) add up the boxes, and then finally (4) divide by the total number of people surveyed. The normal distribution is continuous. Some distributions are *discrete* and don't display this property; there are many examples, but the simplest is maybe the *Poisson* distribution that is only capable of spewing out whole numbers. The problem with statistical distributions is they have a tendency to look terrifyingly complex. I think it's important to remember that each wasn't 'made up' to look a particular way, but rather was mathematically derived to describe samples that were taken under certain conditions. The most important example of this is the normal distribution, which was derived by the 'Prince of Mathematicians' Gauss^[And so the normal distribution is often called the Gaussian distribution; I will use 'normal' because I think it sounds like scary..]. Gauss showed that when you sample something that is affected by lots of independent errors, the distribution of those samples ($x$) will be centred around the true mean ($\mu$^[ 'mu'—the Greek letter for 'm' (mean)]) with a spread defined by the standard deviation ($\sigma$^[ 'sigma'—the Greek letter for 's']) as in this equation: $$ f(x) = \frac{1}{\sqrt{2 \pi} \sigma} e ^\frac{-(x-\mu)^2}{2 \sigma^2} $$ {#eq-normal} This equation defines the normal distribution, and the theorem that distributions will converge on it like this is central to modern statistics and is called the *Central Limit Theorem*. The equation is complicated, and you don't need to learn it, but it's important that you remember that this is why so many histograms tend to look bell-shaped (e.g., @fig-sampling). Things like people's heights are affected by thousands of different factors (*e.g.*, hundreds of genes, dietary choices, lifestyle, and more), and so all those differences add up, through the Central Limit Theorem, to make people's heights normally distributed. Statisticians get so excited about the normal distribution precisely because the Central Limit Theorem points to it: otherwise, it would just be a bizarrely complicated way of drawing a reasonably pretty curve. In @eq-normal above, I referred to a normal distribution's *parameters*: its mean ($\mu$) and standard deviation ($\sigma$). The mean is the centre of the distribution, and most of the samples from that distribution will look more-or-less like that value (it's the distribution's 'central tendency' that, in the 'limit' where you sample from it many times, it will centre on). The standard deviation is the spread of the distribution, and observations from a normal distribution with a greater value will exhibit more variation. You've probably encountered these terms before, but it's critical that you understand there's nothing magical about them. They just happen to be the parameters of the normal distribution; if we were modelling data from another distribution we might use totally different parameters. Indeed, some distributions don't even have defined means or standard deviations! Because frequentists believe that true distributions and probabilities come from sampling, they don't really need maths to find the normal distribution (or any distribution, for that matter). Frequentists believe that, if they keep sampling, the distribution of values they get will converge on the true distribution. Take a look at @fig-sampling, where I show you how much nicer (smoother) and more like the true, underlying distribution, the distribution of samples come when I sample more. If I sampled even more (the grey boxes in @fig-pdf), I could make the width of the histogram bars narrower: I've got more data, so I can split it more finely, and it more accurately shows the shape of the distribution. You would probably agree with me that I could sample even more times, and make the bins even narrower, with the only caveat being I'd have to change the scale on my vertical axis (look at the heights on @fig-sampling; I have to keep changing the limits as I sample more). To derive distributoins, frequentists write a mathematical equation where they sample an infinite number of times. This allows them to make the bars infinitely narrow (because they're doing loads of sampling), and allows them to draw a curve along the mid-points of those bars (they're infinitely narrow, so really they form a curve anyway). So how do we figure out the probability of a particular observation? Remember how height isn't really a single number? How when we thought about people 190 cm tall, we really were asking how many were taller than 185 cm and shorter than 195 cm? It's exacxtly the same with our statistical distributions, becuause to a frequentist the underlying logic and maths is all the same. If we want to know the probability of someone being 190 cm tall, we can read that off our histogram: (1) find the box(es) that cover people between 1.85 and 1.95 meters tall, (2) multiply their heights (the frequency) by their widths to get the total number of people in that range, (3) add up the boxes, and then finally (4) divide by the total number of people surveyed. The only difference now is that the boxes are infinitely narrow, and so it's easier to draw a curve rather than drawing loads of boxes (@fig-pdf). When plotting our curve, just as we had to rescale our histograms, now we have to rescale for infinite sampling, and so divide our frequencies so that the total area under the curve is exactly equal to 1. We need that total area under the curve to be 1 because it's a probability distribution: we need a 100\% probability that something will happen, and everything that can happen needs to be on the curve. We call this plot a _probability density function_ because we use them to get probabilities, just as in histograms we are interested in the area under the curve (that's the 'density' part), and it's defined by an equation we've derived assuming infinite sampling so we can call it a 'curve'. The probability density function for the normal distribution is given to you in @equ-normal. Remember: _a probability is just the area under the curve_. ![**Frequentists derive true statistical distributions by sampling towards infinity.** A frequentist can derive a statistical distribution by deriving, mathematically, what happens when you sample an infinite numbers of times. Our histogram's bins become infinitely narrow, and we can draw a curve across the top of them (the red line at the right). We can calculate the probability of any set of observations (*e.g.*, a height of 190 cm) by measuring the area under the curve between two points that define your question (*e.g.*, the area between 185 and 195 cm).](Figs/01-hist-pdf){#fig-pdf} ## Estimating parameters: maximum likelihood So far we've talked about *drawing* samples from a distribution with known parameters. You may be wondering how, given a set of samples, we can determine the distribution (and parameters) most likely to have generated those samples. We've already talked about estimating the probability of observing a range of possible values (*e.g.*, a height between 185 and 195 cm) and, as outlined above, we can't get a probability of an exact height because it's a continuous variable. We can, however, read off probability density of an exact number (*e.g.*, 190 cm): that's the *likelihood*. For example, plug in an observation ($x$) and parameters ($\mu$ and $\sigma$) to @eq-normal, and voilà, you've got the likelihood that your observation was drawn from that normal distribution. Thus to find out what distribution/parameters a set of samples were drawn from, we find the distribution/parameters that maximise the product of the likelihoods of all the datapoints. We take the product ('product' means 'multiply') because it's the probability we got each value *together*, not independently, and the maximum because bigger means more likely/probably. We call this the *maximum likelihood* solution, or ML estimate of a distribution. Notice that it's an estimate, taken from the sample, of the parameter of the true distribution from which those samples were drawn. The true value of the underlying, true distribution is called the *population* estimate. Estimating all these parameters has a cost. Intuitively, you would probably agree that the more things we estimate, the more underlying data we would need. In statistics, the concept of *degrees of freedom* helps us keep track of this. Imagine a friend sampled ten heights from a normal distribution, told you nine of them, and then offered you a million pounds if you exactly guessed the tenth. You'd say they were being unfair, because there's no way you could get it exactly right—you could get close, but not exactly there. Now imagine that they gave you the mean of all ten numbers as well: you would laugh, because using simple algebra^[Perhaps not simple, but for a million pounds?] you could *exactly* predict the missing value using the formula for a mean^[...$\frac{\sum{x}}{n}$...]. The total degrees of freedom of a random sample is equal to the number of data points in it; once you estimate a parameter from that sample, you 'soak up' (use) a degree of freedom and so you can predict one of those data points from all the others and that estimate. In the previous example, calculating the mean of the same meant there were 9 remaining degrees of freedom ($10 - 1$; 10 data points minus one estimated parameter). Thus you always need at least one data point for each parameter you estimate, and Very Bad Things will happen if you don't have that. My general piece of advice is to make sure you have at least ten data points per estimated parameter, but that's based less on statistical rigour and more on experience. The problem with a single ML solution/estimate is it doesn't give you any idea of the *uncertainty* associated with that estimate. One way to get at this is with *confidence intervals*, often abbreviated to CIs. Confidence intervals allow us to bracket the uncertainty around a parameter estimate that we've found through an ML search. To find them, you have to make a decision: how worried are you about making an *error*? If you want to be really, really certain that you're going to construct some confidence intervals that will contain the true parameter value, then you want to have very wide CIs. The problem is, you won't then have a lot of *power* to know with any certainty what the true value is: the interval will be so wide that it will encompass a lot of possible values. Make sure you take a look at @fig-mle, which gives you an example of what all of this looks like; we'll go through this in more detail in the next chapter. In general, the more data you have, the greater your *power* to detect something that's going on in your data. The more complex your model—the more parameters that you need to build your statistical distribution—the more degrees of freedom you will soak up, and so the more data you will need. A statistical model that misses a real effect (*e.g.*, mood is affected by time of day, but you can't find evidence for it in your data) has made a *Type II error*. Conversely, if you find evidence for an effect that isn't really there (*e.g.*, you suggest there is a link between the MMR vaccine and autism, but really there isn't) you have committed the opposite kind of mistake: a *Type I error*. The way to remember the difference is the story of 'The Boy Who Cried Wolf': first he cries there's a wolf when there isn't one (Type I error), then he cries wolf when there is and no one believes him (they commit a Type II error). ![**A graphical overview of Maximum Likelihood estimation.** As described in the text, it's possible to find the best estimate of a given parameter by calculating the likelihood of the observed data if it were drawn from a distribution parameterised by different estimates. Doing so creates a *likelihood curve*, whose *maximum* represents the best-estimate of that particular parameter (the Maximum Likelihood Estimate, or MLE). It is also possible to construct *confidence intervals* by going a certain distance away from that maximum. Note that, because likelihoods are almost always incredibly small, we work with them on the log scale. Note that it is not coincidence that this likelihood curve looks approximately normal; if you should ever look at an empirical ML curve and it does not, be careful, as it likely indicates something has gone Badly Wrong.](Figs/01-ml.png){#fig-mle} ## The zoo of statistical distributions When I was a graduate student, I used to almost 'collect' statistical distributions, and felt that the more I knew, the 'better' I was at statistics. Sadly, I was wrong^[...as I found out in exams]: it's better to understand a few distributions, and where they come from, than to collect them like imaginary pocket monsters. With that in mind, I'm going to give you a brief tour of some useful statistical distributions, where they come from, and what they can be used for. I'm going to assume that you're comfortable with the concept of a normal distribution, either informally (from the text above) or formally (from the extension section below where we derive it from first principles). I'm also not going to cover the $t$ distribution because it's the basis of the $t$-test, which is what we will be covering extensively in the next chapter. Note also that I'm not giving you nice pictures of what these distributions look like: doing so is one of your exercises below... ### Chi-squared ($\chi^2$) A $\chi^2$ distribution is what you expect to get if you take the variances of normal distributions and sum them up. Let's state that semi-formally^[That phrase 'semi-formally' is doing quite a lot of work here. This is less a definition and more the start of a derivation, to be completely honest, but if you've noticed that then you don't really need to be reading this part of the book anyway...]: $$ Q = \sum \sum (x_k - \bar{x_k})^2 $$ {#eq-chisq} Where Q is a $\chi^2$-distributed variable, $x_k$ is an observation of the $k$th random normally-distributed variable we're working with, and $\bar{x_k}$ is the mean of the $k$th normal distribution. Note that there are two summation signs here because we're doing two sets of summing^[...I'm not putting in the formally required subscripts on the summation signs for simplicity]: once within each random variable, summing the squared deviation from the mean (the variance) for each draw within that sample, and then again to sum across all the variables' squared differences. To calculate the probability density function for this variable, you would carry out that operation across the probability density function of a set of normal distributions. It's a useful distribution to know about because it comes up so much; it's essentially a measure of goodness of fit. Now you know how it's derived from the variances, perhaps that makes some sense, or at least perhaps it will after the next chapter where we will be calculating variances over and over again until you're sick of them. ### F The F-distribution is so-named for Ronald Fisher, who worked extensively on model-testing, and it's incredibly important in model testing and we will spend a lot of time working with it next session. Sadly, Fisher was a horrible person, and notably racist even in a time when racism was rampant and widely accepted, so I don't want to dwell on the distribution's origins anymore. Luckily, the distribution itself is quite simple^[Ahem. Well, it is when you do a simplification along the lines of the trick I used above... This is the last trick in this section, I promise!], albeit looking a bit weird because it is two fractions ($\frac{Q_1}{d_2}$ and $\frac{Q_2}{d_2}$) that are themselves in a larger fraction: $$ F = \frac{\frac{Q_1}{d_1}}{\frac{Q_2}{d_2}} = \frac{Q_1d_2}{Q_2d_1} $$ {#eq-f} Where $Q_1$ and $Q_2$ are two $\chi^2$-distributed variables (see above) and $d_1$ and $d_2$ are the degrees of freedom of those variables. As we'll see next time, this is a very useful distribution because it tells you whether there's anything unusual about the differences in variation of two distributions. It's unusual to see it written out in terms of two sets of fractions, but I encourage you to think of it like that as it helps simplify the concepts we will encounter next time. ### The Poisson distribution The Poisson distribution is the first 'oh wow' distribution for many, because it solves a problem you have been repeatedly told, as an undergraduate, wasn't a problem but in reality actually is a Very Big Deal: nearly-but-not-quite-continuous variables. Not everything in life can take non-integer values: lightning can strike once, twice, thirty times, but not half a time. The Poisson distribution can deal with this, and is predicts only whole numbers: $$ p(x=k) = \frac{\lambda^ke^{-\lambda}}{k!} $$ {#eq-poisson} What this equation means is that the probability that you observe $k$ things (lightning strikes, horse kicks, fish swimming past you) in your random variable $x$ is a function of the natural number $e$ and $\lambda$. $\lambda$^[ 'lambda'—the Greek letter for 'l'] is the only parameter in the Poisson distribution, and is equal to both the mean and the variance of the distribution. We need all the fancy $p(x=k)$ on the left-hand side because we must calculate discrete probabilities for each value of $k$, and there can be no $k$ of 1.5 (only 1 or 2), and thus the Poisson distribution is not really technically a continuous distribution. This all makes the maths of deriving the distribution a bit involved, so we won't dwell on it. ### The Binomial distribution Finally, the workhorse of the bioinformatic world: the Binomial distribution. It's the probability of seeing an event happen a certain number of times given how many times you looked and the underlying probability of it happening. For example, it can give you the probability of getting ten heads when you toss twenty coins. It can give you the probability of seeing 5 diseased individuals in a population of 100 given a background prevlance of 1%. It's your best friend, really it is, and it's so commonly used it's got its own special notation: $$ p(x=k) = \binom{n}{k} p^k (1-p)^{n-k} $$ {#eq-binomial} ... where ... $$ \binom{n}{k} = \frac{n!}{k!(n-k)!} $$ {#eq-binom-meta} ...I know, right? This is definitely the worst one to look at the maths of. Bear with me. It's giving us the probability of detecting $k$ 'successes', where a success is 'the event we're interested in happened' and is quite often, in my experience, something that doesn't feel like a success such as 'someone has a disease'. We are doing the thing (testing a patient, tossing a coin) $n$ times, and our probability of success is $p$. Some people introduce the term $q$ in the above equations to represent the probability of a failure, but I think not writing it out makes the logic of @eq-binomial easier to track through: product of the probability you succeed, the probability you fail (which is just when you don't succeed), and a correction for how many ways there are of doing the two. Let's do it with a worked example: tossing four coins and getting two heads. It's the probability of all the success you would need ($p^k$; two heads are $0.5 \times 0.5 = 0.5^2 = 0.25$) multiplied by the failures you would need ($(1-p)^{n-k}$; two tails are $(1-0.5) \times (1-0.5) = (1-0.5)^2 = 0.5^2 = 0.25$). We then correct for all all the different ways of throwing two heads in a row (where H is 'heads' and T is 'tails': HHTT, HTHT, HTTH, TTHH, THTH, THHT—six ways total): $\frac{n!}{k!(n-k)!} = \frac{4!}{2!(4-2)!} = \frac{4!}{2!2!} = \frac{4 \times 3 \times 2 \times 1}{(2 \times 1) \times (2 \times 1)} = \frac{24}{4} = 6$. In total: $0.25 \times 0.25 \times 6 = 0.375$. Look, it's a lot of maths. But trust me, it's a lot easier than drawing out those blasted branching diagrams you may have been made to do previously at A-level... Oh yes and, of course, will be made to do in the exercises. Sorry about that. ## Extension: Deriving the normal distribution There is no need to read this material, particularly on a first read or if you're not very comfortable with mathematics. If you are, however, you may find this interesting. There are plenty of derivations of the normal distribution online, and if you are interested (and have the mathematical background) I would encourage you to look a couple up. It is possible, however, with the equivalent of an A-level in maths, to understand the important principles of the derivation of the normal distribution and I think it is a useful exercise to go through because it somewhat demystifies what the normal distribution represents. It didn't fall from the sky and it wasn't made up for fun, it follows from relatively simple assumptions. With that in mind, what follows is an attempt to capture the spirit of a full derivation. To say a variable is normally distributed is to say that the frequency of observations drawn from that distribution is lesser the further away from the mean of those observations. Let's write that in maths: $$ \frac{df}{dx} = -v(x - \mu) f(x) $$ {#eq-norm1} Where $f(x)$ is the normal distribution we're trying to derive, $\mu$ is the mean, and $v$^[I'm using non-standard terminology for this constant but you will see why at the end.] is a measure of how fast the frequency falls off away from that mean. The left-hand side of the equation is a derivative: we're talking about the rate of change of the distribution as we move along the values of the observation ($x$), so we're starting from this change (our starting observation that motivated the distribution in the first place) to calculate the actual distribution itself ($f(x)$). There's $x-\mu$ because we are defining the shape of the distribution by how far away we are ($x$) from the mean ($\mu$), and we multiply that by $-v$ because as we 'go away from' the mean (become more positive than the mean) we want the frequency to go down (and go up as we approach it—become less negative). Finally, we multiply everything by $f(x)$ because we want everything to be proportional to the frequency at that point along the distribution: if we have a very low probability, we want to change very slowly because we never want to hit zero and have impossible values, and we want to change fast closer to the centre of the distribution. Equally, if you've been paying attention earlier, you remember that probability distributions need to have an area equal to 1, and so having this in here is convenient in that sense too (but more on that later). OK, I think the above was the most challenging part so if you're still with me well done. Now let's re-arrange a bit to get all our $x$s on one side: $$ \frac{df}{f} = -v(x-\mu) dx $$ {#eq-norm2} That wasn't so bad! We just multiplied by $x$^[Alright, I know calculus doesn't quite work like that, but it really does look that way...]! Now let's get rid of those nasty derivatives by integrating both sides: $$ \int \frac{df}{f} = -v \int (x-\mu) dx $$ {#eq-norm3} The integral on the left-hand side is something you can look up in a table (or derive from first principles if you're that good at maths but that bad at remembering things), and the right-hand side is not so bad if you've done calculus but I do appreciate it's not trivial. Once we've done that we have: $$ \ln f = \frac{-v(x-\mu)^2}{2} + c $$ {#eq-norm4} ...where $c$ is an integration constant as per usual. Let's exponentiate to finish this off and finally write this in terms of $f(x)$: $$ f(x) = c e^{\frac{v}{2}(x-\mu)^2} $$ {#eq-nearly-normal} Great! Now this *nearly* looks like @eq-normal, but it isn't quite there, is it? In particular, we have that nasty integration constant $c$ lying around, and we also have $v$. Let's tackle $c$ first, because it's mathematically the hardest part to deal with (and is where we'll stop with the maths) but is conceptually the simplest. Remember that all probability distributions have to have an area under them equal to 1? This is because *something* has to happen, and so the total realm of possible probabilities has to sum to one. You would probably agree with me that @eq-nearly-normal cannot have a negative value (there is an $e$ in there) unless $c$ is negative, thus $c$ is essentially a normalising constant that determines what the area under that curve will be. Thus the tricky part of the maths of the derivation of the normal distribution is finding a general solution for $c$ such that the integral of @eq-nearly-normal (*i.e.*, the area under the curve) always sums to 1 for any given value of $\mu$ or $v$. That's all there is to it; it happens that the solution is to set $c$ to be $\frac{1}{\sqrt{2\pi\sigma}}$ (I'll get to what $\sigma$ is...) but there's no great magic to it, just a lot of fiddling. So now you know where that normalising constant comes from and what it achieves. So now to $v$. Remember that $v$ is a way of measuring how quickly we swoop up/down when we come closer to or further away from the mean ($\mu$). I suppose, therefore, it a measure of how *variable* the distribution is: larger values of $v$ will mean that the distribution moves more 'quickly' around the mean. Hmmm, so I guess that's a bit silly to call it a measure of how variable the distribution is, then, because bigger values of $v$ will tend to be in distributions that are *less* variable (faster sweeps up/down). Why don't we fix that—it's just a variable, after all—and replace it with its inverse ($\frac{1}{v}$)? That would give us: $$ f(x) = c e^{\frac{1}{2v}(x-\mu)^2} $$ {#eq-norm7} Great! Oh, hang on, 'variable' isn't a good word to use because we're studying random variables and that'll be really confusing. Let's call it *variance* instead. Oh, hang on, we can't use $v$ because we use Greek letters for parameters (like the mean—$\mu$) in statistical distributions because otherwise the maths looks much too approachable and we can't scare graduate students. Also, $v$ looks like the Greek letter $\nu$, which is actually said as an 'n' sound. Oh, stuff it, let's go for something silly—let's call it $\sigma$. Oh wait, hang on, we can never have that variable be negative^[This comes from the maths that I ignore about deriving $c$, but the argument I'm giving here is valid] because otherwise we'll come 'down' towards the mean and then go back 'up', and have a distribution that is weirdly shaped. I guess an easy way to ensure that never happens is just to square it, right, because squared values can never be negative^[Shut up at the back, I know they sometimes can, but we're dealing with real random numbers here, OK? Give me a break...], right? OK, let's call it $\sigma^2$ then: $$ f(x) = c e^{\frac{1}{2\sigma^2}(x-\mu)^2} $$ {#eq-norm8} Oh look! We're done! Hooray! ## Exercises Today's exercises are going to seem a little strange, because we don't know enough practical statistics to start working with real datasets yet. Focus on getting comfortable with the *concepts*, and use these questions as ways to test your understanding. These are the most difficult questions you will face in the course, so don't worry if they seem tough. 1. **Fundamentals of probability definitions.** If you're writing answers to these, do not write pages and pages for these questions: you should be able to answer them in a few sentences. a. Your friend Sarah warns you that there's a 10% chance (in a frequentist sense) that Mount St Helens will erupt tomorrow. Your other friend, Michael, scoffs, and says that her probability makes no sense because the volcano will either erupt tomorrow or it won't, and because there is only one tomorrow there can be no frequentist probability of that. What do you think? b. Your friend Michael approaches you with a problem. He has the average blood pressures of two professors, and wants to statistically test whether the professors' pressures are drawn from a Blaargian or a Floogian distribution. The probabilities each measurement are from a Blaargian distribution are $0.1$ and $0.9$, and the probabilities the pressures are from a Floogian distribution are $0.4$ and $0.6$. Which distribution do you think the data are most likely drawn from? Show your working (hint: if the working isn't trivial, you've gone wrong). c. Michael is beginning to disturb you a little, and wants to keep talking about these blood pressures. He tells you that the Blaargian distribution is defined by one parameter ($\alpha$), while the Floogian distribution is defined by two parameters ($\alpha$ and $\beta$). Does this affect your answer to question 1b above? Why? 2. **Going to the statistical zoo.** The aim of this exercise is not to produce hundreds of plots that you laborious copy-paste into Word, but rather to become familiar with both some of R's basic plotting functions and also the shapes of statistical distributions. Use the code below to generate draws from statistical distributions (which you can lookup in `R` help with something like `?rnorm`), and then sketch the impact of changing different parameters either on a piece of paper or, if it is easier for you, using line drawings in Word or something like that. Make sure you look carefully at the horizontal (the 'x' axis) in these plots, because otherwise you will miss differences between plots. a. The normal distribution has two parameters: $\mu$ and $\sigma$. Use the code below to investigate the impact of changing $\mu$ on the normal distribution, and sketch out your findings. ```{r} # Draw 10,000 times from a normal distribution with # a mean (mu) of 0 and sd (sigma) of 1 # and then make a histogram of those samples hist(rnorm(n=10000, mean=0, sd=1)) # Plot out again, changing the mean hist(rnorm(n=10000, mean=5, sd=1)) # Change the plot so that we can have two plots on the screen # at the same time, then plot the histograms side-by-side par(mfrow=c(1,2)) hist(rnorm(n=10000, mean=0, sd=1), main="mean=0") hist(rnorm(n=10000, mean=5, sd=1), main="mean=5") # Reset the plotting window (clearing the two plots) # so that you can plot one thing at a time again par(mfrow=c(1,1)) ``` b. Use the code above to sketch the impact of changing $\sigma$ (note that `R` uses the standard deviation, $\sigma$, to parameterise the normal and not the variance, $\sigma^2$ but, as one is just the square of the other, it needn't bother you). c. The Poisson distirbution has just one parameter, $\lambda$. It can be plotted out like this: ```{r} hist(rpois(n=10000, lambda=1)) ``` ...sketch out how it changes across $\lambda$. Pay attention to the asymmetry of its distribution. d. The $\chi^2$ distribution is parameterised by the number of distributions it's summing over (there are caveats on this; see the next few chapters). Use the `R` function `rchisq`, and its two arguments `n` and `df`, to explore how the distribution changes across values of `df` such as 1, 5, 10, 50, and 100. Note, again, its asymmetry. Are you beginning to see a pattern in how `R` defines its functions that draw **R**andom numbers from a distribution? e. Drawing random numbers from the Binomial distribution often confuses students because it requires three arguments: the number of times to run the simulation (draw numbers) (`n`), the numbers of times to attempt a success within each draw/run (`size`), and the probability of success (`prob`). Make a sketch of how different values of `size` affect how 'discrete' and 'gappy' the distribution looks, but first follow the outline of how to explore the function given below. ```{r} # Let's toss a coin in the air ten times rbinom(n=10, size=1, prob=0.5) # Let's toss ten coins in the air and count the number of heads rbinom(n=1, size=10, prob=0.5) # Let's toss ten coins in the air and count the number of heads ten times rbinom(n=10, size=10, prob=.5) # ... make sense? Great. If not... Play around a bit, then ask me. ``` f. The F distribution has too many parameters to be easy to sketch, so take five on that one. This question is left deliberately blank. 3. **The Binomial distribution.** You cannot become familiar enough with the Binomial distribution; it really does make your life much easier. While it is always possible to find precise, analytical solutions to Binomial problems, you may find that learning a thing or two about how to simulate draws from it comes in handy from time to time. a. Sarah hands you a bag with five marbles in it, two of which are red and the other three black. She says you can pull a marble out and put it back in the bag (*i.e.*, in statistical parlance you are 'drawing with replacement') three times. She asks to calculate, without using the Binomial distribution but by drawing a probability tree diagram, the probability you will draw three red balls. b. She asks you to now estimate the probability you will draw two red balls, again only using a tree diagram. You are so glad that she helped you with that whole blood pressure situation with Michael that, honestly, you just go along with it. c. Sarah, who makes it clear to you that she finds tree diagrams boring as well, asks you to use the Binomial probability equation your biostatistics lecturer gave you to estimate the answers to (a) and (b) for her. She says her lecturer never gave it to her and she feels hard-done-by, and so would like you to do her a favour. d. After calculating the probabilities by hand, Sarah is not so sure she was that hard-done-by. Sarah shares a trick her lecturer taught her with you, which is to use the Binomial distribution density function (`dbinom`) to estimate the probabilities from parts (a) and (b). She asks if you can do it and, not wanting to look foolish, you glance back at your workings from question 2e and give it a go. e. Sarah's so impressed with your ability to think on your feet, she asks if you can work together to estimate the probability of drawing (with replacement) between 30–50 red balls (inclusive) from a bag with 200 balls in total, of which 30 are red and the remaining 170 are black. She knows lots of ways to do it, but feels there's got to be an easy 'one-liner' in `R` that could do it for her. Hint: Take a look in the help entry for `dbinom` (`?dbinom`). f. Sarah notices a question posted by Michael (who you have both now decided to aviod)written in the margin of her class workbook, wondering whether it matters if half of the black balls in the bag are replaced with blue balls. Would this affect any of your answers? Why (not)? g. Sarah notices a final entry in her workbook from Michael, which appears to be somewhat hastily written and followed by frantic scratchings. It seems to follow a question about whether it would be possible to repeat these exercises 'without replacement': in cases where no ball taken from the bag was returned afterwards. Draw out the probability diagram for part (a) and, on the basis of that, consider whether it was wise of Michael to attempt part (e) using the Binomial distribution.