Now we have covered the core concepts of what a probability is, we can discuss test statistics, using “Student’s” \(t\)-test as an example and then moving on to ANOVAs in general. Test statistics are a fundamental component of statistics, and \(t\)-tests are only the tip of the iceberg. The story of the \(t\)-test is a funny one, and worth bearing in mind in our modern world of closed-access publications and non-disclosure agreements. The test was developed by William Gosset, who worked for the Guinness Brewing Company and had signed a contract stating anything he developed was company property. Unwilling to let what is now one of the major foundations of statistics be hidden away by Guinness, he published under the pseudonym ‘Student’.
Science is all about finding something that we wouldn’t otherwise expect, and the concept of a test statistic was developed to help us measure whether something is surprising to us. Intuitively, you might agree that something was surprising if it reliably differed from our expectations. We can formalise that mathematically as:
\[ t = \frac{o - e}{v} \tag{2.1}\]
Where \(o\) is our observation (e.g., the temperature was \(30^\circ C\)), \(e\) is our expectation (e.g., the temperature yesterday was \(30^\circ C\)), and \(v\) is the overall variation in our observation (e.g., the temperature fluctuated by \(\pm 10^\circ C\) that day). Taking those things together, we can generate a test statistic (\(t\)) that measures how reliable (\(v\); the variation) a difference between our observation and expectation (\(o-e\)) is. In lectures, I make studenst recite “o-e over v” over and over again until they feel sick: learn it, it’s important!
This formalisation of a test statistic underlies essentially all modern frequentist hypothesis-testing. It allows an investigator to set a particular level of \(t\) that they think is sufficiently interesting, and then estimate whether their observation is more or less surprising than that. It is critical that you understand that, to someone using any kind of test statistic, the magnitude of the difference (\(o - e\)) is as important a determinant as the background variation (\(v\)). Thus something that has a very small effect size (\(o-e\); impact on a system) could have a very large test statistic value if the variation was low (\(v\)). Consequently you must never use a test statistic as an estimate of the magnitude of an effect or pattern—there are no exceptions to this, ever1.
So how does one go about deciding whether a particular test statistic value is sufficiently important to warrant further investigation? We compare our test statistic with a distribution of expected values. That distribution is often scaled according to the number of degrees of freedom left over after estimating that value, since we can be more confident of a test statistic that was calculated using more data points (more data, more evidence) and a simpler model (fewer estimated parameters is simpler, and simpler is better—we discuss this below). We can rank that observed test statistic within that distribution: that gives us some degree of certainty as to how unusual a particular value was. It’s much easier to understand this with a concrete example, so let’s move on to our first real test statistic: the \(t\)-test.
A \(t\)-test reveals whether there is evidence that two samples of continuous data are drawn from distributions with different means. The simplest definition of \(t\) is:
\[ t = \frac{\bar{x} - \bar{y}}{SE_{diff}} \tag{2.2}\]
Where \(\bar{x}\) is the mean of one sample (\(x\)) and \(\bar{y}\) is the mean of another (\(y\)). \(\bar{x}-\bar{y}\) is our observation. Our expectation is that the difference between the two means will be nothing (0), and so there is nothing in Equation 2.2 for our expectation because our expectation is nothing! \(SE_{diff}\) is the standard error of the difference in those two means: our measure of how variable the two samples are, and so how much faith we can put in the difference we have observed. For completeness, let’s define it mathematically, but I don’t want you to get stuck on the details:
\[ SE_{diff} = \sqrt{\frac{s^2_x}{n_x} + \frac{s^2_y}{n_y}} \tag{2.3}\]
Where \(s^2_x\) and \(s^2_y\) are the observed variances of \(x\) and \(y\), respectively, and \(n_x\) and \(n_y\) are the number of observations in \(x\) and \(y\), respectively. What matters is that you grasp the intuition of these equations: variation within such samples might be mistaken for variation across samples, and so more variable samples will lead to a greater \(SE_{diff}\) and so a lesser \(t\). Equally, we divide our variances by the number of observations, because if we measure something we’d expect it to vary more—we’re giving it more opportunities for variation.
So we know how to calculate a \(t\)-value; how do we tell whether it’s sufficiently large that we should be surprised by our result? First of all, we have to pick what ‘sufficiently large’ means. Most scientists think that a less than 5% chance of seeing something, by chance, is sufficiently surprising. You may think otherwise, but you must pick before you see the result, otherwise you’re cheating. Because the \(t\)-test is a frequentist test, our result is surprising at the 5% level if it is more extreme than the most extreme 5% of the t-distribution. Remember: in frequentist land, frequency is probability, so there’s only a 1% chance of being more extreme than 1% of the distribution. Another way of phrasing this is that there is a less than 1% chance that we would see a \(t\)-value this extreme if there were no real difference between the two distributions. We call this threshold \(\alpha_{crit}\); if your value were more extreme than 5% of the data, we would say the result is statistically significant at \(\alpha_{5\%}\).
If we were doing this by hand, we would calculate the density of the \(t\)-distribution for a given number of observations and rank our observed \(t\)-value within that distribution (see Figure 2.1 for an example of this). Of course, nowadays we have computers to do this for us, but it’s important you know how this works because you need to know what these numbers mean and they will help you when using ‘bootstrapping’ to generate expected distributions yourself later. Notice that if we were working to \(\alpha_{5\%}\) we would need a \(t\) value that was ranked at or beyond the 2.5th or 97.5th quantile: we would be surprised if the difference were really low or really big, and so we have to split our 5% of ranking between these two expectations. Because we care about both possibilities, this is called a two-tailed test: we are testing both the lower and upper ‘tails’ (sides) of the distribution. If we only cared about one kind of value, we might be able to only check one (perhaps the 5th quantile) and we would be performing a one-tailed test. If you find yourself performing a one-tailed test “to make something significant”, you have almost certainly failed as a scientist and so you should stop what you are doing.
Last time I told you about 95% confidence intervals; these intervals map directly onto \(\alpha_{5\%}\). Thus they come with an important caveat: these CIs map out regions where we would be surprised to see certain differences by chance. But that’s subtly different from giving us a probability that we will detect a true difference if one is present. At \(\alpha_{5\%}\), there is an 80% chance that you will detect a difference between two samples if there is one2. There is also, of course a 5% chance you will make a mistake: 5% of the time you will see a value more extreme than 95% of the distribution, simply by chance!
Let’s start off by simulating some data: two samples drawn from two similar, but slightly different, distributions. We’ll then plot their differences using a boxplot.
exercise <- rnorm(50, mean=2, sd=0.5)
no.exercise <- rnorm(50, mean=1.5, sd=1)
boxplot(exercise, no.exercise, names=c("exercise", "no exercise"))
t.test(exercise, no.exercise)
#>
#> Welch Two Sample t-test
#>
#> data: exercise and no.exercise
#> t = 4.1391, df = 67.848, p-value = 9.848e-05
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#> 0.376049 1.076224
#> sample estimates:
#> mean of x mean of y
#> 1.974944 1.248808The output from our \(t\)-test is quite verbose, but it’s got the most important points that we need from that data: the observed t value, and some degrees of freedom (df) associated with it. We went over degrees of freedom last time; \(t\)-test degrees of freedom, sadly, fall into the category of “quite complicated” because there are multiple kinds of \(t\)-test, each with slightly different ways of dealing with your samples having different variances or sample sizes3. The p-value is probably what you care about most: it’s the ranking of the observed \(t\) value in the \(t\)-distribution, with all the nastiness related to upper and lower tails already taken care of. You also get the means of the two distributions, which is useful.
There is one variant of the \(t\)-test that you will want to get to grips with: the paired t-test. This is test assumes that your samples aren’t independent: perhaps you applied a treatment to subjects and recorded something about them before and after the treatment. If individuals vary naturally, you can use this variation to your advantage, and see if individuals vary within the context of this background variation. Paired \(t\)-tests are, invariably, much more powerful than their non-paired cousins. But you must only use a paired \(t\)-test if your data support it: you can’t use a paired \(t\)-test simply to get a “better” p-value. Let’s start by simulating and then fitting a paired \(t\)-test, and then we’ll see the mathematics are quite straightforward.
# Let's pretend you're happier after class, eh?
base.mood <- rnorm(100)
mood.after.class <- base.mood + .3 + rnorm(100)
# We can't detect a difference normally
t.test(base.mood, mood.after.class)
#>
#> Welch Two Sample t-test
#>
#> data: base.mood and mood.after.class
#> t = -0.66395, df = 175.26, p-value = 0.5076
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#> -0.5050714 0.2507866
#> sample estimates:
#> mean of x mean of y
#> -0.16078376 -0.03364137
# But we can when we control for across-individual variation
t.test(base.mood, mood.after.class, paired=TRUE)
#>
#> Paired t-test
#>
#> data: base.mood and mood.after.class
#> t = -1.1596, df = 99, p-value = 0.249
#> alternative hypothesis: true mean difference is not equal to 0
#> 95 percent confidence interval:
#> -0.3446916 0.0904068
#> sample estimates:
#> mean difference
#> -0.1271424
# ...and this final line shows us what's really going on
t.test(base.mood-mood.after.class)
#>
#> One Sample t-test
#>
#> data: base.mood - mood.after.class
#> t = -1.1596, df = 99, p-value = 0.249
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#> -0.3446916 0.0904068
#> sample estimates:
#> mean of x
#> -0.1271424Bear in mind that these data are random each time, so you might not see such a big effect, but (hopefully) the first test is not significant, while the second and third are. What will definitely be the case each time is that the second and third tests will be identical: that’s because a paired \(t\)-test is just a test of whether the differences between the two distributions are different from zero. More formally, we can define a paired \(t\)-test as:
\[ t_{paired} = \frac{\bar{z}}{\sqrt{\frac{s^2_z}{n_z}}} \]
Where
\[ z = x - y \]
All of the terms make sense if you stare at Equation 2.2 and pretend that, instead of comparing \(x\) with \(y\), you’re comparing \(z\) (\(x-y\); our observation) with 0 (our expectation), and I’ve cascaded through the definition of the standard error (because it’s now simpler).
All the statistical models we’ll deal with from now on have at least two variables: the response variable—the thing we are trying to model—and the explanatory variable(s)—the thing(s) we use to explain the response variable. So, for example, if we wanted to explain the mood of people depending on the day of the year, mood would be our response variable, and would go on the vertical (‘y’4) axis on any plots, and day of year would be the explanatory variable, and would go on the horizontal (‘x’) axis.
In mathematical symbols, a linear regression is the solution to the following equation:
\[ \mu = a + bx \tag{2.4}\]
Where \(\mu\) is the model’s predictions of the response variable, \(a\) is the intercept of the equation, \(b\) is the slope of the equation, and \(x\) is the explanatory variable. Note that, in a paper, you would never write out the solution to this equation in the form \(y = 3.23 \times x + 1.23\), as you may have been taught in high school5. \(a\) and \(b\) are called the model coefficients, parameters, or terms, and our goal in statistics is to find values of these terms that best explain our data. In our first session we discussed how maximum likelihood can be used to find the best values for model coefficients, but that it was important to define precisely what we meant by ‘best’. In this case, ‘best’ means values that minimise the squared error of our model, which as defined as:
\[ SSE = \sum(y - \mu)^2 \tag{2.5}\]
Where \(SSE\) is the sum of the squared error of our model, \(\sum\) means “add up (sum) all of the”, \(y\) is the response variable, and \(\mu\) is defined above. This equation is very important, and you should memorise it: it means that a linear regression minimises the error in our model, which is the same as the sum of the (squared) lengths of the red lines in Figure 2.2. The model’s error for each prediction is often called the residuals, and is defined formally as \(y - \mu\). Notice that we try and minimise the sum of the squared errors: this is because otherwise, when summing things up, negative residuals could make it look as if we’d made no error—we need to square the errors to ensure that both over- and under-predicting is equally bad6.
ANOVA stands for ANalysis Of VAriance, and differs from linear regression only in that the explanatory variables are discrete (often called factors in R). Thus an ANOVA could be used to test whether mood varies among two different groups, such as people who exercise (group a) and people who don’t (group b). The definition of SSE, which is still the quantity we are minimising, remains the same, but the definition of our predictions changes somewhat:
\[ \mu = \bar{y_{group}} \tag{2.6}\]
Where \(\mu\) is our predicted response, as before, and \(\bar{y_{group}}\) is the mean value of our response variable in a group7. Considering it graphically is often the easiest way to think about it, so let’s consider it with reference to Figure 2.2 b. Our response variable is split into groups according to the explanatory variable, and in this case there are two groups. Our ANOVA is trying to find an estimate for each group that minimises the difference between the values of the response in those groups and the group-level estimate—it’s minimising the sums of the squared errors, exactly like in regression. Handily, it just so happens that the estimate that minimises the squared error will always be the mean value for the groups. Actually, this is also the case for regression (because, remember, they’re the same thing): a regression line will always pass through the mean value of the response and explanatory variables.
It might seem confusing that ANOVA is called ANOVA, given that we use ANOVA to estimate whether treatments vary in their mean values. The reason is that all the significance values (see below) are based around variance: it’s the variation, or variance, that our model does or doesn’t explain, in comparison with the overall variance in our response variable, that tells us whether we’ve found a statistically significant difference or not.
You now know how regressions and ANOVAs are fit to data, and where all the estimates come from. What you don’t know is how to test model significance. As we discussed at the beginning of the course, model significance is a tricky subject: you don’t want to get too obsessed with whether something is or is not past an arbitrary threshold. Regardless, there is value to ascribing how ‘surprising’ your results are, and this is a very meaningful concept that a significance test can help you get at.
All data varies: ideally a statistical model would explain lots of that variation, and thus the left-over residual variation would be quite small. If we’re going to know whether our model is doing a good job, we should compare the amount of variation it explains with the overall amount of variation in our data, right? In ANOVA, we define the overall variation in our data in the following way:
\[ SST = \sum (y - \bar{y})^2 \]
Where \(SST\) is the total variation, \(y\) is (still) our response variable, and \(\bar{y}\) is the mean of that response variable. This is ‘just another’ sum of squares (hence SS-T), and Figure 2.3 shows what it looks like in the contexts of the regression and ANOVA you’ve already encountered. Let’s write an equation to formalise the relationship between variance attributed to our model (\(SSM\)8), variance in our total data (\(SST\)), and our residual variation (\(SSE\); defined above):
\[ SST = SSM + SSE \tag{2.7}\]
This makes sense, right? All the variation in our data can be partitioned into that explained by our model (\(SSM\)), and that that isn’t (\(SSE\)). The variation our model explains is the total variation in the data minus whatever left-over variation the model doesn’t explain. It also means, quite handily, that we can calculate how much variation is explained by our model using \(SSE\) and \(SST\), both of which are quite easy to calculate. It’s a simple re-arrangement of Equation 2.7 above, but here it is again just for completeness:
\[ SSM = SST - SSE \]
So now we’re ready to, finally, figure out whether our model is interesting. I’d define ‘interesting’ as “explains a lot of variation compared with what it doesn’t explain, given how many parameters are in my model”. I’ve added that bit about parameters in there because, as we discussed in our first lecture, a really complicated model should do a better job of explaining things, but it isn’t necessarily parsimonious. Let’s write that out in maths:
\[ F = \frac{SSM \times (n-k-1)}{SSE \times (k-1)} \tag{2.8}\]
Where \(n\) is the number of data points, and \(k\) is the number of parameters in our model (other than the overall mean, \(\bar{y}\), which is used to calculate \(SST\)). \(F\) is a number that represents that ratio between the variances, and, just like the \(t\)-statistic, we have null expectations of how big we would expect it to be by chance for a dataset of a given size (\(n\)) with a given number of parameters (\(k\)). This \(F\) statistic is used to generate our \(p\)-value, and so determine whether our model is interesting enough to be ‘statistically significant’. Just as the \(t\) distribution we compared our observed \(t\) value to was defined by the degrees of freedom left once \(t\) has been calculated, so too our \(F\) distribution is defined by two parameters: the degrees of freedom on the top and bottom of its equation [(\(n-k-1\)) and (\(k-1\))]. You might like to quickly look back at the previous session and you will see that this definition of an F-statistic lines up with the general definition you were given earlier.
What’s wonderful about this framework is that you can use it to compare any set of models. When we put \(SSE\) on the bottom of the \(F\) calculation, we’re comparing our model (the alternative model or alternative hypothesis) to a null model (or null hypothesis) where there is no difference between the treatments. If we wanted to, we could pick a different model, calculate its \(SSM\) (using the equations above), and then compare those two models. This way you can take a very complicated model and ask if it explains significantly more variation than would be expected by chance than a simpler model. This kind of nested testing is a useful approach, and it’s something we’ll be coming back to in later lectures.
Phew! That’s a lot of theory, so let’s put it into practice. Remember that we will be returning to all of this over the coming sessions, so if it seems a bit overwhelming now just keep calm and carry on. You’ll have plenty of time to get to grips with all of this. Let’s start off by simulating some data, and then see how we would analyse that in a regression. We’re going to simulate 100 data-points according to a regression equation with an intercept (\(a\)) of 2 and a slope (\(b\)) of 1.5.
cont.explanatory <- rnorm(100)
response <- 2 + 1.5*cont.explanatory + rnorm(100, sd=0.5)
plot(response ~ cont.explanatory)
reg.model <- lm(response ~ cont.explanatory)
summary(reg.model)
#>
#> Call:
#> lm(formula = response ~ cont.explanatory)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.64614 -0.32838 0.02891 0.37422 0.95669
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.02371 0.05342 37.88 <2e-16 ***
#> cont.explanatory 1.42392 0.05165 27.57 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.5341 on 98 degrees of freedom
#> Multiple R-squared: 0.8858, Adjusted R-squared: 0.8846
#> F-statistic: 760.1 on 1 and 98 DF, p-value: < 2.2e-16There’s a lot of output in the summary table here, so let’s focus only on the ones that we know about right now (eventually, you will understand everything here). The (Intercept) estimate is our value of \(a\), and the cont.explanatory estimate is our value of \(b\)—they should be close to 2 and 1.5. Don’t worry about the p-values and standard errors just yet—we’ll get to them next time. R summaries our residuals for us, so that we can check that they are normally distributed and centred at zero (more on another time, but they don’t look like they aren’t). At the very bottom is the F-statistic: we get a value for it (it’s very big, because the model explains a lot of variation compared with what it doesn’t) and we get some degrees of freedom which map onto the \(n\) and \(k\) values we saw above: we estimated 1 extra parameter (the slope), and because there are 100 data-points and we estimated two parameters in total (the intercept and the slope) that is listed as \(98\).
An ANOVA is exactly the same maths as a regression, only our explanatory variables are discrete, not continuous. So let’s split our data into three groups according to the continuous explanatory variable cont.explanatory and fit an ANOVA to that data instead.
disc.explanatory <- ifelse(cont.explanatory < -1, "a",
ifelse(cont.explanatory < 1, "b", "c")
)
anova.model <- lm(response ~ disc.explanatory)
summary(anova.model)
#>
#> Call:
#> lm(formula = response ~ disc.explanatory)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.3253 -0.4147 0.0392 0.5349 2.1164
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.09588 0.18381 -0.522 0.603
#> disc.explanatoryb 2.28521 0.21268 10.745 <2e-16 ***
#> disc.explanatoryc 4.28851 0.27482 15.605 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.8423 on 97 degrees of freedom
#> Multiple R-squared: 0.7188, Adjusted R-squared: 0.713
#> F-statistic: 124 on 2 and 97 DF, p-value: < 2.2e-16Everything’s the same as it was before except for how the estimates for the levels within the factor disc.explanatory. The first ‘level’ within our your factor—a—is the elephant in the room: it’s not named, but it’s there. Remember that all of our models fit an intercept, and that’s why it doesn’t count as an extra parameter9? The first treatment level is reported as the (Intercept), and R’s estimate of its mean is roughly 0. disc.explanatoryb is not the mean for level b: it is the difference between a and b. It is a contrast, so-called because it contrasts the effect of a and b. Similarly, disc.explanatoryc is also a contrast, describing the difference between a and c. Never forget that these are not estimates of the mean, but rather estimates of the contrasts.
If you want to look ‘behind the curtain’, you can use the function summary.aov to get the summary of your models in terms of the sums of squares (as defined above). What’s confusing in R is that it gives you different summaries by default depending on whether you used the lm function to fit a regression/ANOVA or the aov function. You can switch between them at will, as I show below. Remember: it’s easy to do this because a regression and an ANOVA are the same thing.
model <- lm(response ~ disc.explanatory)
summary(model)
#>
#> Call:
#> lm(formula = response ~ disc.explanatory)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.3253 -0.4147 0.0392 0.5349 2.1164
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.09588 0.18381 -0.522 0.603
#> disc.explanatoryb 2.28521 0.21268 10.745 <2e-16 ***
#> disc.explanatoryc 4.28851 0.27482 15.605 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.8423 on 97 degrees of freedom
#> Multiple R-squared: 0.7188, Adjusted R-squared: 0.713
#> F-statistic: 124 on 2 and 97 DF, p-value: < 2.2e-16
summary.lm(model) # identical to above
#>
#> Call:
#> lm(formula = response ~ disc.explanatory)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.3253 -0.4147 0.0392 0.5349 2.1164
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.09588 0.18381 -0.522 0.603
#> disc.explanatoryb 2.28521 0.21268 10.745 <2e-16 ***
#> disc.explanatoryc 4.28851 0.27482 15.605 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.8423 on 97 degrees of freedom
#> Multiple R-squared: 0.7188, Adjusted R-squared: 0.713
#> F-statistic: 124 on 2 and 97 DF, p-value: < 2.2e-16
summary.aov(model) # different
#> Df Sum Sq Mean Sq F value Pr(>F)
#> disc.explanatory 2 175.95 87.97 124 <2e-16 ***
#> Residuals 97 68.82 0.71
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model <- aov(response ~ disc.explanatory)
summary.aov(model) # identical to above
#> Df Sum Sq Mean Sq F value Pr(>F)
#> disc.explanatory 2 175.95 87.97 124 <2e-16 ***
#> Residuals 97 68.82 0.71
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.lm(anova.model) # different
#>
#> Call:
#> lm(formula = response ~ disc.explanatory)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.3253 -0.4147 0.0392 0.5349 2.1164
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -0.09588 0.18381 -0.522 0.603
#> disc.explanatoryb 2.28521 0.21268 10.745 <2e-16 ***
#> disc.explanatoryc 4.28851 0.27482 15.605 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.8423 on 97 degrees of freedom
#> Multiple R-squared: 0.7188, Adjusted R-squared: 0.713
#> F-statistic: 124 on 2 and 97 DF, p-value: < 2.2e-16Finally, how can we compare different models to see which is the more parsimonious, or explains the data better? By using the anova function, which performs comparisons between different models using the logic I describe above. We’re going to go through this in much more detail in the next chapter, so focus on the pick-picture for the time-being and don’t worry about the details just ye. Let’s use anova to compare between our two different model formulations, and then test an overly complicated model to see what extra it gives us (not very much).
# Test our two initial models
anova(reg.model, anova.model)
#> Analysis of Variance Table
#>
#> Model 1: response ~ cont.explanatory
#> Model 2: response ~ disc.explanatory
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 98 27.954
#> 2 97 68.824 1 -40.87
# Make some random data
silly.explanatory <- rnorm(100)
# Fit a silly model
silly.model <- lm(response ~ silly.explanatory + cont.explanatory)
# Test whether our silly model does better than our other model
anova(silly.model, reg.model)
#> Analysis of Variance Table
#>
#> Model 1: response ~ silly.explanatory + cont.explanatory
#> Model 2: response ~ cont.explanatory
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 97 27.793
#> 2 98 27.954 -1 -0.161 0.5619 0.4553When we perform an ANOVA of our more complicated model (two slope parameters—we’ll worry about how you can get multiple slopes into a regression next time) versus our simpler model, we find a non-significant (greater than 5%) p-value. This means there’s no evidence to choose one model over the other on the basis of the variance they explain, given how many parameters are in each model. So which model do we think is the more interesting? The one with fewer parameters, because it’s the simplest.
There’s more to statistical modelling than finding the best model of a set (the relative performance of your model); you need to know how well your model performs overall (its absolute performance). The classic way to calculate this is with something called an \(r^2\), which measures how much of the overall variation in your data your model explains. This provides a measure of the degree of fit of our model to our data: a better fit means it explains the data better. We can define \(r^2\) as:
\[ r^2 = \frac{SSM}{SST} \tag{2.9}\]
Which is, frankly, just a fancy way of saying that the \(r^2\) is a number between 1 and 0: it’s 0 if your model explains no variation in your data, and 1 if your model explains all the variation in your data. Thus good models will have higher \(r^2\) values. The astute of you will have noticed that I’ve still not told you what a good model is, and the truth is that there is no answer: different fields have different criteria. In macro-ecology, an \(r^2\) of 15% (0.15) is often considered a pretty good model; in epidemiology, an \(r^2\) of 15% could get someone killed.
A model with more parameters should do a better job of fitting the data simply by chance, and so many scientists prefer the concept of an adjusted \(r^2\). The most common kind, which is by default calculated by R, is defined as:
\[ r^2_{adjusted} = \frac{SSM \times (n - 1)}{SST \times (n - k - 1)} \tag{2.10}\]
Where \(k\) is the number of parameters in your model, and \(n\) is the number of data points in your dataset. This number penalises models with more parameters, and makes it possible to (sort-of) directly compare models. Personally, I’m not a fan. The correction is intended to be like the correction in an ANOVA (compare the brackets in Equation 2.10 with those in Equation 2.8); there’s not much more logic to it than this, and if you want to calculate something that’s like an ANOVA… why not just fit an ANOVA? I think you’re better off just bearing in mind how complicated your model is when looking at the \(r^2\) value, and for anything else framing a hypothesis test along the lines I show you in the next chapter. Regardless of how you feel, everyone will assume you’re reporting a vanilla \(r^2\) value unless you give them reason to otherwise, so make sure you’re clear if you report the \(r^2_{adjusted}\) in a manuscript. A thing to note is that, in all R model summaries, the \(r^2\) as I’ve defined it above is termed the multiple \(r^2\) in order to distinguish it from the adjusted \(r^2\). Do not confuse the two!
Model criticism means different things to different people, but all of them have in common a desire to make sure that your model is fitting your data well, and that predictions from that model are reasonable. Thus \(r^2\) and the deviance statistics we’ll discuss next time are valid forms of model criticism: they tell you whether your model is explaining a lot of the variation in your model.
Statisticians are often additionally concerned about whether the assumptions underlying their model have been met. For example, a standard linear model assumes that its residuals are normally distributed: it models \(y \sim normal(\mu, \sigma^2)\), and so the error about our mean prediction (\(\mu\)) should be normally distributed with a constant variance (\(\sigma^2\)). Generalised Linear Models (and mixed effects models, which we’ll come to next time) make different assumptions, and it’s important to think about what those assumptions are and whether your model appears to be obeying them.
Let’s start off by looking at what a well-fitting linear model “looks like”. The code below simulates some data and generates the plots you can see in Figure 2.4; it makes use of a function called plt.mod to make all four model diagnostic plots appear in the same window. You can display each in turn by running something like plot(model); try it for yourself if you’re curious. Don’t focus on the details of the R code, though, as I don’t expect you to understand what’s going on with par and friends; use the function as you like, but focus on the statistics, not the programming.
# Function to make plots look prettier
plt.mod <- function(model, ...){
old.par <- par(mfrow=c(2,2))
plot(model, ...)
par(old.par)
}
# Simulate some well-behaved data
x <- rnorm(100)
y <- x + rnorm(100, sd=.5)
# Fit model and plot its diagnostic plots
model <- lm(y ~ x)
plt.mod(model)
Figure 2.4 shows what the diagnostic plots for this model, whose assumptions have been met, look like. The legend of Figure 2.4 goes through what these mean in detail, and so you should read through that legend carefully now. Now I’m going to show you how to simulate some data that do not meet the assumptions of a standard linear regression, and show you what the diagnostic plots look like for both an lm and glm fit to those data.
# Simulate some poorly-behaved data
x <- rnorm(100)
y <- numeric(100)
for(i in seq_along(x))
y[i] <- rpois(1, abs(x[i]))
# Fit a bad model and plot its diagnostic plots
bad.model <- lm(y ~ x)
plt.mod(bad.model)
Figure 2.5 shows the model diagnostic plots for the bad.model—a model whose assumptions do not make it a good fit for those data. Again, the legend for Figure 2.5 goes through what each plot shows in some detail, and so you should read through that legend carefully before returning to the main text.
Linear models assume that the residuals—the left-over variation that our model doesn’t predict—are normally distributed and have constant variance. Looking at Figure 2.4 and Figure 2.5, you can see that the top- and bottom-left plots reveal the non-constant variance in the residuals. In particular, the bottom-left plot shows a very strong upward trend in the data with non-constant variance. Perhaps most-importantly, however, it’s clear that the residuals aren’t even normally distributed: while the top-right plot in Figure 2.4 shows a nice match of the residuals to what we would expect of them if they were normally distributed, in Figure 2.5 the extremes of the residuals are vastly too large. Indeed, the discontinuities suggest the underlying data aren’t even continuous and so neither are the residuals (they aren’t, they’re drawn from the Poisson distribution). Finally, the bottom-right plot of Figure 2.4 shows no trend in the outliers (points with undue influence on our model), whereas Figure 2.5 shows that some over-predictions are biasing our results. This is in keeping with our residuals being too extreme in the model: these outliers are so off that they’re pulling our model away from the ‘true’ results.
Model criticism is hard, and it’s definitely true to say that your get better at it with practice. This is something I would really encourage you to go over with your supervisor if you’re feeling confused as they will definitely have seen their fair share of residual plots. While there are statistical tests you can do to see if your model’s residuals are normally distributed10, I’d advise you not to spend much time on them. The more tests you run, the greater your chance of making a mistake because of the multiple testing and \(\alpha_{crit}\) problems we’ve already discussed. Moreover, fundamentally, you need to be sure that the assumptions of your test of your test are met. How do you know if they are? Well, fit another test. But wait, how do you know if that test of the test’s test’s assumptions were met? Well, fit another test. But wait, how do you know if… Fundamentally, at some level all statistics requires a modicum of subjectivity: in Bayesian statistics it gets pulled in earlier in the form of prior belief (see session one), and you should count yourself lucky that you’re only encountering it now.
So, what should you do if your diagnostic plots don’t look good? First off, try fitting a more complex model by adding more terms into the model to explain what’s going on. Second, if that doesn’t work, try a more complex model statistically, such as a GLM or something even more complex (see our next session). Remember that, if fitting a GLM, the diagnostic plots in R will look the same for that GLM as they do here: the onus will be on you to know what the assumptions of that new test are, and to update your expectations accordingly. I cover both of these options in the coming sessions. If you want my advice, though, I would say this: it is better to add biological reality than statistical complexity. Add more variables to your model, and make those variables mean something biologically. Even the simplest of statistical models, with sufficient biological information added, can reveal something useful about your system.
As ever, below is some code to load today’s dataset into R, and then the questions themselves.
# Load the mileage data like this
data <- read.csv("auto-mpg-simplified.csv", as.is=TRUE)
# Plot the data out like this
with(data, plot(mpg ~ acceleration))
with(data, boxplot(mpg ~ cylinders))
# Load the crab data in like this
data <- read.csv("crabs_simplified.csv", as.is=TRUE)
# Plot the data out like this
with(data, boxplot(diameter ~ sex, main="Are male crabs larger?"))
with(data, boxplot(diameter ~ source, main="Are tourist-caught crabs larger?"))
us in a specialised crab trap12 or by a tourist (them) on with a lure on the end of a piece of fishing line. Each crab has a sex (male/female) and a measured carapace width.continuous.variable ~ grouping.variable…in much the same way that I show you how to use boxplot above. c. Use a \(t\)-test to check whether we were catching crabs that were of a different size to those of the tourists. Why do you think this might (not) be? d. The local tourist industry of Wells-next-the-Sea relies on tourists coming year-on-year to catch crabs. What do you think the results you’ve found might mean for the sustainability of the economy of Wells-next-the-Sea? Note: I’m not looking for a great long diatribe, rather a few short, snappy sentences considering what might happen if tourists catch all the large crabs… 3. The following questions require you to use a (slightly simplified) version of data from Quinlan (1993; Proceedings on the Tenth International Conference of Machine Learning, 236–243). This dataset describes the fuel efficiency of several kinds of car (mpg—miles per gallon), and can be described as a function of various other properties of the cars. a. Using the code above, plot a scatter-plot of each car’s fuel efficiency as a function of its acceleration. Do you think there is a relationship between the two? b. Fit a regression to the data to test whether there is a correlation between fuel efficiency and acceleration rate. c. Fit an ANOVA to the data to test whether cars with an unusually large number of cylinders have a higher fuel efficiency. Do you think this relationship is correlation or causation? d. Compare the results of your ANOVA to that of a \(t\)-test. Do they give the same answer13? e. The ‘average’ car emits \(8,887\) grams of CO\(_2\) per gallon burned14. Estimate how many miles you drive per year, and, using the mean of all cars in this dataset (mean function in R), calculate how many kilograms of CO\(_2\) you would release in the ‘average’ car every year. Compare the emissions from the average fuel consumption with those of a car with a fuel efficiency 10mpg higher than this average. f. The ‘average’ price for a gallon of gas is \(\$2.54\). How much money would you spend on gas per year in the ‘average’ car? How much would you spend in our hypothetical car with its slight fuel efficiency increase? 4. To help you master your statistical concepts at the same time as your programming, below are three fill-in-the-blanks R functions that test your understanding of the definitions of sums of squares and test statistics. I would strongly encourage you to see if you can fill in these blanks using the information in the handout above. You don’t need to necessarily understand how to program to fill in these blanks—look at the equations and think what you might fill in.
my.t.test <- function(x, y){
observation <- mean(___) - mean(___)
std.err <- sqrt((sd(___)^2)/length(___) + (sd(___)^2)/length(___))
t <- ___ / ___
return(t)
}
calc.sst <- function(x){
____ <- mean(x)
____ <- (____ - x)^2
return(____)
}
calc.ssm <- function(resp, explan){
model <- lm(___ ~ ___)
sst <- calc.sst(___)
residuals <- residuals(model)
sse <- sum((____)^2)
ssm <- ___ - ___
return(ssm)
}
r2 <- function(resp, explan){
sst <- ____
sse <- ____
____ <- ___ / ___
return(____)
}Maximum likelihood by hand. There is nothing magical about how maximum likelihood (ML15) estimates are derived. They are found by a computing trying, in careful order, different parameter values in ways that are only slightly mathematically cleverer than what you are about to do now. The purpose of this exercise is to demystify that entire process for you.
Michael is acting oddly again, and has now measured the average systolic blood pressures of 5 professors (he assures you they were taken ‘at rest’ and for reasons you can’t quite put your finger on that concerns you): 120, 130, 110, 105, and 121. Using the function dnorm, estimate the likelihood that these numbers were drawn from a normal distribution with a mean of 120 and standard deviation of 5. The code below may help, but remember that dnorm will return the density (likelihood) of a given set of input data and so, to calculate the overall likelihood of all the data, you will have to do something to all of those likelihoods. Hint: likelihood is just a fancy word for probability, and so how do we PROD probabilities to get their joint probability16?
dnorm(x=c(-5, 5), mean=0, sd=2)
#> [1] 0.00876415 0.00876415Michael is unsatisfied, and decides he will have to keep taking measurements until he can be sure. You convince him that a 95% confidence interval will tell him all he needs to know about how variable the blood pressure measurements could be. His eyes light up, and so you decide to fix the estimate of the standard deviation of the distribution at the standard deviation of your data (sd(c(120,130,110,105,121))) and use that to find (roughly) the 95% CI. Find (to the nearest whole number) the upper and lower limits of the mean. Hint: look above at how I define the 95% CIs in log-likelihood units, and remember that the R function for the natural logarithm is log.
Michael is becoming agitated, and says that your estimate is wrong because you assumed the estimate of the SD was the same as the sample’s SD. Reassure him by now fixing your estimate of the mean to be what you found in part (b), and estimate the 95% CIs on the standard deviation.
Michael’s eyes light up, as he realises you have essentially implemented the first pass through something called Brent’s method (you don’t care, you just want him out of your office). He asks you one final thing before he leaves: what is the probability, given all this, that he would find a person with a systolic blood pressure below 100? You remember that you can find this for, say, a value of -1 from a ‘standard normal’ distribution (mean of 0 and standard deviation of 1) with pnorm(-1, mean=0, sd=1) but can you do it for Michael’s value17?
Michael leaves, but you still feel uneasy. You check your own systolic blood pressure and find it’s 150. Given you know, from (d), how to find the area to the left of a distribution, can you figure out how to estimate the probability that you would have a resting blood pressure of at least 150? As a bonus exercise, use this information to reflect upon whether you should be spending as much time with Michael as you do, and whether you should tell someone about him.
For an empirical example of this problem see Pearse et al. (2013; Ecology, 94(12): 2861.↩︎
The derivation of this 80% power requires some knowledge of linear algebra, and so I’m not going to go through it. If you are interested, the book Likelihood in the further reading is a good introduction to this. You should also note that 80% power is the expectation under most conditions; many tests have a greater power than this, many also have much less, so buyer beware… In the next chapter’s exercises, I give you an introduction to how to estimate power for any statistical test.↩︎
You really don’t want to know the details, I assure you. Because this is a fundamental aspect of statistics, the implications of a tiny change here are quite large in other areas, so the whole field can look quite navel-gaze-y unless you know a lot more statistics than you do right now.↩︎
Please, please, please: try not to use the ‘y’ and ‘x’ axis terminology. Some people are very snobby about it, and those people can get quite annoyed about it. I’m not annoyed by it, but sometimes appearances matter.↩︎
There are examples of how to present statistical results at the end of your handout.↩︎
OK, there’s a bit more going on here… Students are often confused as to why we couldn’t ‘simply’ use the absolute error, i.e., minimise \(|x - \mu|\). The answer is that we can, but we were only able to figure out mathematically what that would mean in the last few decades. The reason we do it is not that by doing so we happen to ensure normally-distributed residuals; squaring something doesn’t guarantee that the output is normally distributed (indeed, arguably the opposite because normal distributions have no upper or lower bounds and if you square something the result is bounded at zero). The real reason comes from the derivation of the normal distribution you’ve already seen: the maths was worked through and, if the residuals are normally distributed, then minimising the sum of the squared errors will give you the solution to Equation 2.5. The squaring is because, in Equation 1.1, such a squared term is found: there’s a fundamental link between the mathematical equations.↩︎
If this seems like a strange equation, that’s because it is. If you knew linear algebra, I could write the equation for \(\mu\) in exactly the same way for a regression and an ANOVA (remember: they’re the same thing), but as you don’t I can’t. Mathematicians have this awful habit of writing very short equations and acting as if the fact they’re short makes them simple; sorry!↩︎
Beware: different textbooks use different names for the things I’ve described above. So, for example, in “The R Book”, \(SSM\) is called \(SSA\). The names don’t matter, but the concepts do. If you’re happy with the concepts, don’t worry about anything else.↩︎
You can force R not to do this, and fit a model without an intercept. Running lm(response ~ disc.explanatory -1) would do that. By all means use it to help clear up confusions for yourself but you must never report or trust estimates of confidence or significance from such a model. Once you understand what -1 does, and have read to the end of the handout, see if you can figure out what anova(reg.model, lm(response ~ 1)) tests.↩︎
shapiro.test↩︎
Pearse et al. 2014; arXiv:1404.0290.↩︎
A glorified bucket↩︎
If you’re interested, stare for a few moments at the equation for a \(t\)-test and ask yourself if you can see any fundamental differences between it and an ANOVA. Ask me if you would like to ponder this a while…↩︎
https://www.epa.gov/greenvehicles/greenhouse-gas-emissions-typical-passenger-vehicle.↩︎
Later in the course, when we cover machine learning, you are going to see just how awfully confusing it is that ML as an acronym is used for both maxmimum likelihood and machine learning. Not my fault, don’t shoot the messenger.↩︎
Even bigger hint: to get the probability that two things happen, you multiply them together. To get the probability that three will happen, you multiply them together. Another word for multiplying is taking the ‘product’ of something. Does this PROD you into thinking about an R function you could use to help with this? PROD PROD? PROD?…↩︎
Are you seeing that pattern in the naming ofr R probability functions yet?…↩︎