12 Correlation and linear regression
The goal in this chapter is to introduce correlation and linear regression. These are the standard tools that statisticians rely on when analysing the relationship between continuous predictors and continuous outcomes.
12.1 Correlations
In this section we’ll talk about how to describe the relationships between variables in the data. To do that, we want to talk mostly about the correlation between variables. But first, we need some data (Table 12.1).
12.1.1 The data
Let’s turn to a topic close to every parent’s heart: sleep. The data set we’ll use is fictitious, but based on real events. Suppose I’m curious to find out how much my infant son’s sleeping habits affect my mood. Let’s say that I can rate my grumpiness very precisely, on a scale from 0 (not at all grumpy) to
Next, I’ll take a look at some basic descriptive statistics and, to give a graphical depiction of what each of the three interesting variables looks like, Figure 12.1 plots histograms. One thing to note: just because jamovi can calculate dozens of different statistics doesn’t mean you should report all of them. If I were writing this up for a report, I’d probably pick out those statistics that are of most interest to me (and to my readership), and then put them into a nice, simple table like the one in Table 12.1.2 Notice that when I put it into a table, I gave everything “human readable” names. This is always good practice. Notice also that I’m not getting enough sleep. This isn’t good practice, but other parents tell me that it’s pretty standard.
12.1.2 The strength and direction of a relationship
We can draw scatterplots to give us a general sense of how closely related two variables are. Ideally though, we might want to say a bit more about it than that. For instance, let’s compare the relationship between baby.sleep and dani.grump (Figure 12.2 (a)), left, with that between dani.sleep and dani.grump (Figure 12.2 (b)), right. When looking at these two plots side by side, it’s clear that the relationship is qualitatively the same in both cases: more sleep equals less grump! However, it’s also pretty obvious that the relationship between dani.sleep and dani.grump is stronger than the relationship between baby.sleep and dani.grump. The plot on the right is “neater” than the one on the left. What it feels like is that if you want to predict what my mood is, it’d help you a little bit to know how many hours my son slept, but it’d be more helpful to know how many hours I slept.
In contrast, let’s consider the two scatterplots shown in Figure 12.3. If we compare the scatterplot of “baby.sleep v dani.grump” (left) to the scatterplot of “baby.sleep v dani.sleep” (right), the overall strength of the relationship is the same, but the direction is different. That is, if my son sleeps more, I get more sleep (positive relationship, right-hand side), but if he sleeps more then I get less grumpy (negative relationship, left-hand side).
12.1.3 The correlation coefficient
We can make these ideas a bit more explicit by introducing the idea of a correlation coefficient (or, more specifically, Pearson’s correlation coefficient), which is traditionally denoted as r. The correlation coefficient between two variables
[Additional technical detail3]
By standardising the covariance, not only do we keep all of the nice properties of the covariance discussed earlier, but the actual values of r are on a meaningful scale: r = 1 implies a perfect positive relationship and
12.1.4 Calculating correlations in jamovi
Calculating correlations in jamovi can be done by clicking on the ‘Regression’ – ‘Correlation Matrix’ button. Transfer all four continuous variables across into the box on the right to get the output in Figure 12.5.
12.1.5 Interpreting a correlation
Naturally, in real life you don’t see many correlations of
Correlation | Strength | Direction |
---|---|---|
-1.00 to -0.90 | Very strong | Negative |
-0.90 to -0.70 | Strong | Negative |
-0.70 to -0.40 | Moderate | Negative |
-0.40 to -0.20 | Weak | Negative |
-0.20 to 0.00 | Negligible | Negative |
0.00 to 0.20 | Negligible | Positive |
0.20 to 0.40 | Weak | Positive |
0.40 to 0.70 | Moderate | Positive |
0.70 to 0.90 | Strong | Positive |
0.90 to 1.00 | Very strong | Positive |
Note that I say a rough guide. There aren't hard and fast rules for what counts as strong or weak relationships. It depends on the context |
However, something that can never be stressed enough is that you should always look at the scatterplot before attaching any interpretation to the data. A correlation might not mean what you think it means. The classic illustration of this is “Anscombe’s Quartet” (Anscombe, 1973), a collection of four data sets. Each data set has two variables, an
You’d think that these four data sets would look pretty similar to one another. They do not. If we draw scatterplots of
12.1.6 Spearman’s rank correlations
The Pearson correlation coefficient is pretty useful, but it does have shortcomings. One issue stands out: what it actually measures is the strength of the linear relationship between two variables. In other words, what it gives you is a measure of the extent to which the data all tend to fall on a single, perfectly straight line. Often, this is a pretty good approximation to what we mean when we say “relationship”, and so the Pearson correlation is a good thing to calculate. Sometimes though, it isn’t.
One very common situation where the Pearson correlation isn’t quite the right thing to use arises when an increase in one variable
To illustrate, consider the data plotted in Figure 12.7, showing the relationship between hours worked and grade received for 10 students taking some class. The curious thing about this (highly fictitious) data set is that increasing your effort always increases your grade. It might be by a lot or it might be by a little, but increasing effort will never decrease your grade. If we run a standard Pearson correlation, it shows a strong relationship between hours worked and grade received, with a correlation coefficient of
How should we address this? Actually, it’s really easy. If we’re looking for ordinal relationships all we have to do is treat the data as if it were ordinal scale! So, instead of measuring effort in terms of “hours worked”, lets rank all
Okay, so let’s have a look at our students when we rank them from worst to best in terms of effort and reward Table 12.3.
Hmm. These are identical. The student who put in the most effort got the best grade, the student with the least effort got the worst grade, etc. As the table above shows, these two rankings are identical, so if we now correlate them we get a perfect relationship, with a correlation of 1.0.
What we’ve just re-invented is Spearman’s rank order correlation, usually denoted
12.2 Scatterplots
Scatterplots are a simple but effective tool for visualising the relationship between two variables, like we saw with the figures in the section on Correlations. It’s this latter application that we usually have in mind when we use the term “scatterplot”. In this kind of plot each observation corresponds to one dot. The horizontal location of the dot plots the value of the observation on one variable, and the vertical location displays its value on the other variable. In many situations you don’t really have a clear opinion about what the causal relationship is (e.g., does A cause B, or does B cause A, or does some other variable C control both A and B). If that’s the case, it doesn’t really matter which variable you plot on the x-axis and which one you plot on the y-axis. However, in many situations you do have a pretty strong idea which variable you think is most likely to be causal, or at least you have some suspicions in that direction. If so, then it’s conventional to plot the cause variable on the x-axis, and the effect variable on the y-axis. With that in mind, let’s look at how to draw scatterplots in jamovi, using the same parenthooddata set (i.e. parenthood.csv) that I used when introducing correlations.
Suppose my goal is to draw a scatterplot displaying the relationship between the amount of sleep that I get (dani.sleep) and how grumpy I am the next day (dani.grump). There are two different ways in which we can use jamovi to get the plot that we’re after. The first way is to use the ‘Plot’ option under the ‘Regression’ – ‘Correlation Matrix’ button, giving us the output shown in Figure 12.8. Note that jamovi draws a line through the points, we’ll come onto this a bit later in the section on What is a linear regression model?. Plotting a scatterplot in this way also allows you to specify ‘Densities for variables’ and this option adds a density curve showing how the data in each variable is distributed.
The second way do to it is to use one of the jamovi add-on modules. This module is called ‘scatr’ and you can install it by clicking on the large ‘
12.2.1 More elaborate options
Often you will want to look at the relationships between several variables at once, using a scatterplot matrix (in jamovi via the ‘Correlation Matrix’ – ‘Plot’ command). Just add another variable, for example baby.sleep to the list of variables to be correlated, and jamovi will create a scatterplot matrix for you, just like the one in Figure 12.10.
12.3 What is a linear regression model?
Stripped to its bare essentials, linear regression models are basically a slightly fancier version of the Pearson correlation (see Correlations), but they are actually much more powerful tools. We’ll return to the parenthood.csv file that we were using to illustrate how correlations work. Recall that, in this data set we were trying to find out why Dani is so very grumpy all the time and our working hypothesis was that I’m not getting enough sleep. We drew a scatterplots to help us examine the relationship between the amount of sleep I get and my grumpiness the following day, as in Figure 12.9, and as we saw that this corresponded to a correlation of
This is not highly surprising. The line that I’ve drawn in Figure 12.11(b) doesn’t “fit” the data very well, so it doesn’t make a lot of sense to propose it as a way of summarising the data, right? This is a very simple observation to make, but it turns out to be very powerful when we start trying to wrap just a little bit of maths around it. To do so, let’s start with a refresher of some high school maths. The formula for a straight line is usually written like this:
The two variables are
Hmm. Looks like the same formula, but there’s some extra frilly bits in this version. Let’s make sure we understand them. Firstly, notice that I’ve written
Excellent, excellent. Next, I can’t help but notice that, regardless of whether we’re talking about the good regression line or the bad one, the data don’t fall perfectly on the line. Or, to say it another way, the data
12.4 Estimating a linear regression model
Okay, now let’s redraw our pictures but this time I’ll add some lines to show the size of the residual for all observations. When the regression line is good, our residuals (the lengths of the solid black lines) all look pretty small, as shown in Figure 12.12(a), but when the regression line is a bad one the residuals are a lot larger, as you can see from looking at Figure 12.12(b). Hmm. Maybe what we “want” in a regression model is small residuals. Yes, that does seem to make sense. In fact, I think I’ll go so far as to say that the “best fitting” regression line is the one that has the smallest residuals. Or, better yet, since statisticians seem to like to take squares of everything why not say that:
The estimated regression coefficients,
and , are those that minimise the sum of the squared residuals, which we could either write as or as .
Yes, yes that sounds even better. And since I’ve indented it like that, it probably means that this is the right answer. And since this is the right answer, it’s probably worth making a note of the fact that our regression coefficients are estimates (we’re trying to guess the parameters that describe a population!), which is why I’ve added the little hats, so that we get
At this point, we now have a concrete definition for what counts as our “best” choice of regression coefficients,
12.4.1 Linear regression in jamovi
To run my linear regression, open up the ‘Regression’ – ‘Linear Regression’ analysis in jamovi, using the parenthood.csv data file. Then specify dani.grump as the ‘Dependent Variable’ and dani.sleep as the variable entered in the ‘Covariates’ box. This gives the results shown in Figure 12.13, showing an intercept
12.4.2 Interpreting the estimated model
The most important thing to be able to understand is how to interpret these coefficients. Let’s start with
12.5 Multiple linear regression
The simple linear regression model that we’ve discussed up to this point assumes that there’s a single predictor variable that you’re interested in, in this case dani.sleep. In fact, up to this point every statistical tool that we’ve talked about has assumed that your analysis uses one predictor variable and one outcome variable. However, in many (perhaps most) research projects you actually have multiple predictors that you want to examine. If so, it would be nice to be able to extend the linear regression framework to be able to include multiple predictors. Perhaps some kind of multiple regression model would be in order?
Multiple regression is conceptually very simple. All we do is add more terms to our regression equation. Let’s suppose that we’ve got two variables that we’re interested in; perhaps we want to use both dani.sleep and baby.sleep to predict the dani.grump variable. As before, we let
As before,
12.5.1 Doing it in jamovi
Multiple regression in jamovi is no different to simple regression. All we have to do is add additional variables to the ‘Covariates’ box in jamovi. For example, if we want to use both dani.sleep and baby.sleep as predictors in our attempt to explain why I’m so grumpy, then move baby.sleep across into the ‘Covariates’ box alongside dani.sleep. By default, jamovi assumes that the model should include an intercept. The coefficients we get this time are shown in Table 12.4.
The coefficient associated with dani.sleep is quite large, suggesting that every hour of sleep I lose makes me a lot grumpier. However, the coefficient for baby.sleep is very small, suggesting that it doesn’t really matter how much sleep my son gets. What matters as far as my grumpiness goes is how much sleep I get. To get a sense of what this multiple regression model looks like, Figure 12.14 shows a 3D plot that plots all three variables, along with the regression model itself.
[Additional technical detail7]
12.6 Quantifying the fit of the regression model
So we now know how to estimate the coefficients of a linear regression model. The problem is, we don’t yet know if this regression model is any good. For example, the regression.1 model claims that every hour of sleep will improve my mood by quite a lot, but it might just be rubbish. Remember, the regression model only produces a prediction
12.6.1 The value
Once again, let’s wrap a little bit of mathematics around this. Firstly, we’ve got the sum of the squared residuals:
which we would hope to be pretty small. Specifically, what we’d like is for it to be very small in comparison to the total variability in the outcome variable:
While we’re here, let’s calculate these values ourselves, not by hand though. Let’s use something like Excel or another standard spreadsheet programme. I have done this by opening up the parenthood.csv file in Excel and saving it as parenthood rsquared.xls so that I can work on it. The first thing to do is calculate the
- Create a new column called ’ Y.pred ’ using the formula ’ = 125.97 + (-8.94
dani.sleep) ’. - Calculate the
by creating a new column called ’ (Y-Y.pred)^2 ’ using the formula ’ = (dani.grump - Y.pred)^2 ’. - Then, at the bottom of this column calculate the sum of these values, i.e. ’ sum( ( Y-Y.pred)^2 ) ’.
- At the bottom of the dani.grump column, calculate the mean value for dani.grump (NB Excel uses the word ’ AVERAGE ’ rather than ’ mean ’ in its function).
- Then create a new column, called ’ (Y - mean(Y))^2 ) ’ using the formula ’ = (dani.grump - AVERAGE(dani.grump))^2 ’.
- Then, at the bottom of this column calculate the sum of these values, i.e. ‘sum( (Y - mean(Y))^2 )’.
- Calculate
by typing into a blank cell the following: ‘= 1 - (SS(resid) / SS(tot) )’.
This gives a value for
Naturally, you don’t actually need to type all these commands into Excel yourself if you want to obtain the
12.6.2 The relationship between regression and correlation
At this point we can revisit my earlier claim that regression, in this very simple form that I’ve discussed so far, is basically the same thing as a correlation. Previously, we used the symbol
12.6.3 The adjusted value
One final thing to point out before moving on. It’s quite common for people to report a slightly different measure of model performance, known as “adjusted
[Additional technical detail9]
This adjustment is an attempt to take the degrees of freedom into account. The big advantage of the adjusted
An obvious question then is whether you should report
12.7 Hypothesis tests for regression models
So far we’ve talked about what a regression model is, how the coefficients of a regression model are estimated, and how we quantify the performance of the model (the last of these, incidentally, is basically our measure of effect size). The next thing we need to talk about is hypothesis tests. There are two different (but related) kinds of hypothesis tests that we need to talk about: those in which we test whether the regression model as a whole is performing significantly better than a null model, and those in which we test whether a particular regression coefficient is significantly different from zero.
12.7.1 Testing the model as a whole
Okay, suppose you’ve estimated your regression model. The first hypothesis test you might try is the null hypothesis that there is no relationship between the predictors and the outcome, and the alternative hypothesis that the data are distributed in exactly the way that the regression model predicts.
[Additional technical detail10]
We’ll see much more of the
12.7.2 Tests for individual coefficients
The
I can’t help but notice that the estimated regression coefficient for the baby.sleep variable is tiny (
How can we test this? Well, if the central limit theorem is kind to us we might be able to guess that the sampling distribution of
I’ll skip over the reasons why, but our degrees of freedom in this case are
In any case, this
12.7.3 Running the hypothesis tests in jamovi
To compute all of the statistics that we have talked about so far, all you need to do is make sure the relevant options are checked in jamovi and then run the regression. If we do that, as in Figure 12.15, we get a whole bunch of useful output.
The ‘Model Coefficients’ at the bottom of the jamovi analysis results shown in Figure 12.15 provides the coefficients of the regression model. Each row in this table refers to one of the coefficients in the regression model. The first row is the intercept term, and the later ones look at each of the predictors. The columns give you all of the relevant information. The first column is the actual estimate of
The only thing that the coefficients table itself doesn’t list is the degrees of freedom used in the
12.8 Regarding regression coefficients
Before moving on to discuss the assumptions underlying linear regression and what you can do to check if they’re being met, there’s two more topics I want to briefly discuss, both of which relate to the regression coefficients. The first thing to talk about is calculating confidence intervals for the coefficients. After that, I’ll discuss the somewhat murky question of how to determine which predictor is most important.
12.8.1 Confidence intervals for the coefficients
Like any population parameter, the regression coefficients
[Additional technical detail13]
In jamovi we had already specified the ‘95% Confidence interval’ as shown in Figure 12.15, although we could easily have chosen another value, say a ‘99% Confidence interval’ if that is what we decided on.
12.8.2 Calculating standardised regression coefficients
One more thing that you might want to do is to calculate “standardised” regression coefficients, often denoted
The basic idea is quite simple; the standardised coefficients are the coefficients that you would have obtained if you’d converted all the variables to z-scores before running the regression.14 The idea here is that, by converting all the predictors to z-scores, they all go into the regression on the same scale, thereby removing the problem of having variables on different scales. Regardless of what the original variables were, a
[Additional technical detail15]
To make things even simpler, jamovi has an option that computes the
These results clearly show that the dani.sleep variable has a much stronger effect than the baby.sleep variable. However, this is a perfect example of a situation where it would probably make sense to use the original coefficients b rather than the standardised coefficients
12.9 Assumptions of regression
The linear regression model that I’ve been discussing relies on several assumptions. In Model checking we’ll talk a lot more about how to check that these assumptions are being met, but first let’s have a look at each of them.
- Linearity. A pretty fundamental assumption of the linear regression model is that the relationship between
and actually is linear! Regardless of whether it’s a simple regression or a multiple regression, we assume that the relationships involved are linear. - Independence: residuals are independent of each other. This is really just a “catch all” assumption, to the effect that “there’s nothing else funny going on in the residuals”. If there is something weird (e.g., the residuals all depend heavily on some other unmeasured variable) going on, it might screw things up. Independence isn’t something that you can check directly and specifically with diagnostic tools, but if your regression diagnostics are messed up then think carefully about the independence of your observations and residuals.
- Normality. Like many of the models in statistics, basic simple or multiple linear regression relies on an assumption of normality. Specifically, it assumes that the residuals are normally distributed. It’s actually okay if the predictors
and the outcome variables are non-normal, so long as the residuals are normal. See the Checking the normality of the residuals section. - Equality (or “homogeneity”) of variance. Strictly speaking, the regression model assumes that each residual
is generated from a normal distribution with mean 0, and (more importantly for the current purposes) with a standard deviation that is the same for every single residual. In practice, it’s impossible to test the assumption that every residual is identically distributed. Instead, what we care about is that the standard deviation of the residual is the same for all values of , and (if we’re being especially diligent) all values of every predictor in the model.
So, we have four main assumptions for linear regression (that neatly form the acronym LINE). And there are also a couple of other things we should also check for:
- Uncorrelated predictors. The idea here is that, in a multiple regression model, you don’t want your predictors to be too strongly correlated with each other. This isn’t “technically” an assumption of the regression model, but in practice it’s required. Predictors that are too strongly correlated with each other (referred to as “collinearity”) can cause problems when evaluating the model. See Checking for collinearity section.
- No “bad” outliers. Again, not actually a technical assumption of the model (or rather, it’s sort of implied by all the others), but there is an implicit assumption that your regression model isn’t being too strongly influenced by one or two anomalous data points because this raises questions about the adequacy of the model and the trustworthiness of the data in some cases. See the section on Outliers and anomalous data.
12.10 Model checking
The main focus of this section is regression diagnostics, a term that refers to the art of checking that the assumptions of your regression model have been met, figuring out how to fix the model if the assumptions are violated, and generally to check that nothing “funny” is going on. I refer to this as the “art” of model checking with good reason. It’s not easy, and while there are a lot of easily available tools that you can use to diagnose and maybe even cure the problems that affect your model (if there are any, that is!), you really do need to exercise a certain amount of judgement when doing this.
In this section I describe several different things you can do to check that your regression model is doing what it’s supposed to. It doesn’t cover the full space of things you could do, but it’s still much more detailed than what is often done in practice – unfortunately! But it’s important that you get a sense of what tools are at your disposal, so I’ll try to introduce a bunch of them here. Finally, I should note that this section draws quite heavily from Fox & Weisberg (2011), the book associated with the “car” package that is used to conduct regression analysis in
12.10.1 Three kinds of residuals
The majority of regression diagnostics revolve around looking at the residuals, and there are several different kinds of residual that we might consider. In particular, the following three kinds of residuals are referred to in this section: “ordinary residuals”, “standardised residuals”, and “Studentised residuals”. There is a fourth kind that you’ll see referred to in some of the Figures, and that’s the “Pearson residual”. However, for the models that we’re talking about in this chapter the Pearson residual is identical to the ordinary residual.
The first and simplest kind of residuals that we care about are ordinary residuals. These are the actual raw residuals that I’ve been talking about throughout this chapter so far. The ordinary residual is just the difference between the predicted value
This is of course what we saw earlier, and unless I specifically refer to some other kind of residual, this is the one I’m talking about. So there’s nothing new here. I just wanted to repeat myself. One drawback to using ordinary residuals is that they’re always on a different scale, depending on what the outcome variable is and how good the regression model is. That is, unless you’ve decided to run a regression model without an intercept term, the ordinary residuals will have mean 0 but the variance is different for every regression. In a lot of contexts, especially where you’re only interested in the pattern of the residuals and not their actual values, it’s convenient to estimate the standardised residuals, which are normalised in such a way as to have a standard deviation of 1.
[Additional technical detail16]
The third kind of residuals are Studentised residuals (also called “jackknifed residuals”) and they’re even fancier than standardised residuals. Again, the idea is to take the ordinary residual and divide it by some quantity in order to estimate some standardised notion of the residual.17
Before moving on, I should point out that you don’t often need to obtain these residuals yourself, even though they are at the heart of almost all regression diagnostics. Most of the time the various options that provide the diagnostics, or assumption checks, will take care of these calculations for you. Even so, it’s always nice to know how to actually get hold of these things yourself in case you ever need to do something non-standard.
12.10.2 Checking the linearity of the relationship
We should check for the linearity of the relationships between the predictors and the outcomes. There’s a few different things that you might want to do in order to check this. Firstly, it never hurts to just plot the relationship between the predicted values
In any case, in order to get a more detailed picture it’s often more informative to look at the relationship between the predicted values and the residuals themselves. Again, in jamovi you can save the residuals to the data set and then draw a scatterplot of the predicted values against the residual values, as in Figure 12.18. As you can see, not only does it draw the scatterplot showing the predicted value against the residuals, you can also plot a line through the data that shows the relationship between the two. Ideally, this should be a straight, perfectly horizontal line. In practice, we’re looking for a reasonably straight or flat line. This is a matter of judgement.
More advanced versions of the same plot are produced by checking ‘Residuals plots’ in the regression analysis ‘Assumption checks’ options in jamovi. These are useful for checking linearity, normality and equality of variance assumptions, and we look at these in more detail in Section 12.10.3. This option not only draws plots comparing the predicted values to the residuals, it does so for each individual predictor.
12.10.3 Checking the normality of the residuals
Like many of the statistical tools we’ve discussed in this book, regression models rely on a normality assumption. In this case, we assume that the residuals are normally distributed. The first thing we can do is draw a QQ-plot via the ‘Assumption Checks’ – ‘Assumption Checks’ – ‘Q-Q plot of residuals’ option. The output is shown in Figure 12.19, showing the standardised residuals plotted as a function of their theoretical quantiles according to the regression model.
Another thing we should check is the relationship between the predicted (fitted) values and the residuals themselves. We can get jamovi to do this using the ‘Residuals Plots’ option, which provides a scatterplot for each predictor variable, the outcome variable, and the predicted values against residuals, see Figure 12.20. In these plots we are looking for a fairly uniform distribution of dots, with no clear bunching or patterning of the dots. Looking at these plots, there is nothing particularly worrying as the dots are fairly evenly spread across the whole plot. There may be a little bit of non-uniformity in plot (b), but it is not a strong deviation and probably not worth worrying about.
If we were worried, then in a lot of cases the solution to this problem (and many others) is to transform one or more of the variables. We discussed the basics of variable transformation in Section 6.3, but I do want to make special note of one additional possibility that I didn’t explain fully earlier: the Box-Cox transform. The Box-Cox function is a fairly simple one and it’s very widely used.18
You can calculate it using the BOXCOX function in the ‘Compute’ variables screen in jamovi.
12.10.4 Checking equality of variance
The regression models that we’ve talked about all make an equality (i.e.homogeneity) of variance assumption: the variance of the residuals is assumed to be constant. To plot this in jamovi first we need to calculate the square root of the (absolute) size of the residual19 and then plot this against the predicted values, as in Figure 12.21. Note that this plot actually uses the standardised residuals rather than the raw ones, but it’s immaterial from our point of view. What we’re looking to see here is a straight, horizontal line running through the middle of the plot.20
12.10.5 Checking for collinearity
Another regression diagnostic is provided by variance inflation factors (VIFs), which are useful for determining whether or not the predictors in your regression model are too highly correlated with each other. There is a variance inflation factor associated with each predictor
If you’ve only got two predictors, the VIF values are always going to be the same, as we can see if we click on the ‘Collinearity’ checkbox in the ‘Regression’ – ‘Assumptions’ options in jamovi. For both dani.sleep and baby.sleep the VIF is
To give a sense of how we could end up with a model that has bigger collinearity problems, suppose I were to run a much less interesting regression model, in which I tried to predict the day on which the data were collected, as a function of all the other variables in the data set. To see why this would be a bit of a problem, let’s have a look at the correlation matrix for all four variables (Figure 12.22).
We have some fairly large correlations between some of our predictor variables! When we run the regression model and look at the VIF values, we see that the collinearity is causing a lot of uncertainty about the coefficients. First, run the regression, as in Figure 12.23 and you can see from the VIF values that, yep, that’s some mighty fine collinearity there.
12.10.6 Outliers and anomalous data
One danger that you can run into with linear regression models is that your analysis might be disproportionately sensitive to a smallish number of “unusual” or “anomalous” observations. I discussed this idea previously in Section 5.2.3 in the context of discussing the outliers that get automatically identified by the boxplot option under ‘Exploration’ – ‘Descriptives’, but this time we need to be much more precise. In the context of linear regression, there are three conceptually distinct ways in which an observation might be called “anomalous”. All three are interesting, but they have rather different implications for your analysis.
The first kind of unusual observation is an outlier. The definition of an outlier (in this context) is an observation that is very different from what the regression model predicts. An example is shown in Figure 12.24, the outlier has an unusual value on the outcome (y-axis location), but not the predictor (x-axis location) and lies a long way from the regression line. In practice, we operationalise this concept by saying that an outlier is an observation that has a very large residual,
Outliers are interesting: a big outlier might correspond to junk data, e.g., the variables might have been recorded incorrectly in the data set, or some other defect may be detectable. Note that you shouldn’t throw an observation away just because it’s an outlier. But the fact that it’s an outlier is often a cue to look more closely at that case and try to find out why it’s so different.
The second way in which an observation can be unusual is if it has high leverage, which happens when the observation is very different from all the other observations. This doesn’t necessarily have to correspond to a large residual. If the observation happens to be unusual on all variables in precisely the same way, it can actually lie very close to the regression line. An example of this is shown in Figure 12.25. The leverage of an observation is operationalised in terms of its hat value, usually written
In general, if an observation lies far away from the other ones in terms of the predictor variables, it will have a large hat value (as a rough guide, high leverage is when the hat value is more than 2-3 times the average; and note that the sum of the hat values is constrained to be equal to
This brings us to our third measure of unusualness, the influence of an observation. A high influence observation is an outlier that has high leverage. That is, it is an observation that is very different to all the other ones in some respect, and also lies a long way from the regression line. This is illustrated in Figure 12.26. Notice the contrast to the previous two figures. Outliers don’t move the regression line much and neither do high leverage points. But something that is both an outlier and has high leverage, well that has a big effect on the regression line. That’s why we call these points high influence, and it’s why they’re the biggest worry. We operationalise influence in terms of a measure known as Cook’s distance.23
In order to have a large Cook’s distance an observation must be a fairly substantial outlier and have high leverage. As a rough guide, Cook’s distance greater than 1 is often considered large (that’s what I typically use as a quick and dirty rule).
In jamovi, information about Cook’s distance can be calculated by clicking on the ‘Cook’s Distance’ checkbox in the ‘Assumption Checks’ – ‘Data Summary’ options. When you do this, for the multiple regression model we have been using as an example in this chapter, you get the results as shown in Figure 12.27.
You can see that, in this example, the mean Cook’s distance value is
An obvious question to ask next is, if you do have large values of Cook’s distance what should you do? As always, there’s no hard and fast rule. Probably the first thing to do is to try running the regression with the outlier with the greatest Cook’s distance24 excluded and see what happens to the model performance and to the regression coefficients. If they really are substantially different, it’s time to start digging into your data set and your notes that you no doubt were scribbling as your ran your study. Try to figure out why the point is so different. If you start to become convinced that this one data point is badly distorting your results then you might consider excluding it, but that’s less than ideal unless you have a solid explanation for why this particular case is qualitatively different from the others and therefore deserves to be handled separately.
12.11 Model selection
One fairly major problem that remains is the problem of “model selection”. That is, if we have a data set that contains several variables, which ones should we include as predictors, and which ones should we not include? In other words, we have a problem of variable selection. In general, model selection is a complex business but it’s made somewhat simpler if we restrict ourselves to the problem of choosing a subset of the variables that ought to be included in the model. Nevertheless, I’m not going to try covering even this reduced topic in a lot of detail. Instead, I’ll talk about two broad principles that you need to think about, and then discuss one concrete tool that jamovi provides to help you select a subset of variables to include in your model. First, the two principles:
It’s nice to have an actual substantive basis for your choices. That is, in a lot of situations you the researcher have good reasons to pick out a smallish number of possible regression models that are of theoretical interest. These models will have a sensible interpretation in the context of your field. Never discount the importance of this. Statistics serves the scientific process, not the other way around.
To the extent that your choices rely on statistical inference, there is a trade off between simplicity and goodness of fit. As you add more predictors to the model you make it more complex. Each predictor adds a new free parameter (i.e., a new regression coefficient), and each new parameter increases the model’s capacity to “absorb” random variations. So the goodness of fit (e.g.,
) continues to rise, sometimes trivially or by chance, as you add more predictors no matter what. If you want your model to be able to generalise well to new observations you need to avoid throwing in too many variables.
This latter principle is often referred to as Ockham’s razor and is often summarised in terms of the following pithy saying: do not multiply entities beyond necessity. In this context, it means don’t chuck in a bunch of largely irrelevant predictors just to boost your
In any case, what we need is an actual mathematical criterion that will implement the qualitative principle behind Ockham’s razor in the context of selecting a regression model. As it turns out there are several possibilities. The one that I’ll talk about is the Akaike information criterion (AIC) (Akaike, 1974) simply because it’s available as an option in jamovi.25
The smaller the AIC value, the better the model performance. If we ignore the low level details it’s fairly obvious what the AIC does. On the left we have a term that increases as the model predictions get worse; on the right we have a term that increases as the model complexity increases. The best model is the one that fits the data well (low residuals, left-hand side) using as few predictors as possible (low K, right-hand side). In short, this is a simple implementation of Ockham’s razor.
AIC can be added to the ‘Model Fit Measures’ output Table when the ‘AIC’ checkbox is clicked, and a rather clunky way of assessing different models is seeing if the AIC value is lower if you remove one or more of the predictors in the regression model. This is the only way currently implemented in jamovi, but there are alternatives in other more powerful programmes, such as R. These alternative methods can automate the process of selectively removing (or adding) predictor variables to find the best AIC. Although these methods are not implemented in jamovi, I will mention them briefly below just so you know about them.
12.11.1 Backward elimination
In backward elimination you start with the complete regression model, including all possible predictors. Then, at each “step” we try all possible ways of removing one of the variables, and whichever of these is best (in terms of lowest AIC value) is accepted. This becomes our new regression model, and we then try all possible deletions from the new model, again choosing the option with lowest AIC. This process continues until we end up with a model that has a lower AIC value than any of the other possible models that you could produce by deleting one of its predictors.
12.11.2 Forward selection
As an alternative, you can also try forward selection. This time around we start with the smallest possible model as our start point, and only consider the possible additions to the model. However, there’s one complication. You also need to specify what the largest possible model you’re willing to entertain is.
Although backward and forward selection can lead to the same conclusion, they don’t always.
12.11.3 A caveat
Automated variable selection methods are seductive things, especially when they’re bundled up in (fairly) simple functions in powerful statistical programmes. They provide an element of objectivity to your model selection, and that’s kind of nice. Unfortunately, they’re sometimes used as an excuse for thoughtlessness. No longer do you have to think carefully about which predictors to add to the model and what the theoretical basis for their inclusion might be. Everything is solved by the magic of AIC. And if we start throwing around phrases like Ockham’s razor, well it sounds like everything is wrapped up in a nice neat little package that no-one can argue with.
Or, perhaps not. Firstly, there’s very little agreement on what counts as an appropriate model selection criterion. When I was taught backward elimination as an undergraduate, we used
What does this mean in practice? Well, you could go and spend several years teaching yourself the theory of model selection, learning all the ins and outs of it so that you could finally decide on what you personally think the right thing to do is. Speaking as someone who actually did that, I wouldn’t recommend it. You’ll probably come out the other side even more confused than when you started. A better strategy is to show a bit of common sense. If you’re staring at the results of an automated backwards or forwards selection procedure, and the model that makes sense is close to having the smallest AIC but is narrowly defeated by a model that doesn’t make any sense, then trust your instincts. Statistical model selection is an inexact tool, and as I said at the beginning, interpretability matters.
12.11.4 Comparing two regression models
An alternative to using automated model selection procedures is for the researcher to explicitly select two or more regression models to compare to each other. You can do this in a few different ways, depending on what research question you’re trying to answer. Suppose we want to know whether or not the amount of sleep that my son got has any relationship to my grumpiness, over and above what we might expect from the amount of sleep that I got. We also want to make sure that the day on which we took the measurement has no influence on the relationship. That is, we’re interested in the relationship between baby.sleep and dani.grump, and from that perspective dani.sleep and day are nuisance variable or covariates that we want to control for. In this situation, what we would like to know is whether dani.grump ~ dani.sleep + day + baby .sleep (which I’ll call Model 2, or M2) is a better regression model for these data than dani.grump ~ dani.sleep + day (which I’ll call Model 1, or M1). There are two different ways we can compare these two models, one based on a model selection criterion like AIC, and the other based on an explicit hypothesis test. I’ll show you the AIC based approach first because it’s simpler, and follows naturally from discussion in the last section. The first thing I need to do is actually run the two regressions, note the AIC for each one, and then select the model with the smaller AIC value as it is judged to be the better model for these data. Actually, don’t do this just yet. Read on because there is an easy way in jamovi to get the AIC values for different models included in one table.26
A somewhat different approach to the problem comes out of the hypothesis testing framework. Suppose you have two regression models, where one of them (Model 1) contains a subset of the predictors from the other one (Model 2). That is, Model 2 contains all of the predictors included in Model 1, plus one or more additional predictors. When this happens we say that Model 1 is nested within Model 2, or possibly that Model 1 is a submodel of Model 2. Regardless of the terminology, what this means is that we can think of Model 1 as a null hypothesis and Model 2 as an alternative hypothesis. And in fact we can construct an
Okay, so that’s the hypothesis test that we use to compare two regression models to one another. Now, how do we do it in jamovi? The answer is to use the ‘Model Builder’ option and specify the Model 1 predictors dani.sleep and day in ‘Block 1’ and then add the additional predictor from Model 2 (baby.sleep) in ‘Block 2’, as in Figure 12.28. This shows, in the ‘Model Comparisons’ Table, that for the comparisons between Model 1 and Model 2,
We can also use this ‘Model Comparison’ option to display a table that shows the AIC and BIC for each model, making it easy to compare and identify which model has the lowest value, as in Figure 12.28.
12.12 Summary
- Want to know how strong the relationship is between two variables? Calculate Correlations.
- Drawing Scatterplots.
- Basic ideas about What is a linear regression model? and Estimating a linear regression model.
- Multiple linear regression.
- Quantifying the fit of the regression model using
. - Hypothesis tests for regression models.
- In Regarding regression coefficients we talked about calculating Confidence intervals for the coefficients and Calculating standardised regression coefficients.
- The Assumptions of regression and Model checking.
- Regression Model selection.
I’ve noticed that in jamovi you can also specify an ‘ID’ variable type, but for our purposes it does not matter how we specify the ID variable as we won’t be including it in any analyses.↩︎
Actually, even that table is more than I’d bother with. In practice most people pick one measure of central tendency, and one measure of variability only.↩︎
The formula for the Pearson’s correlation coefficient can be written in several different ways. I think the simplest way to write down the formula is to break it into two steps. Firstly, let’s introduce the idea of a covariance. The covariance between two variables
and is a generalisation of the notion of the variance and is a mathematically simple way of describing the relationship between two variables that isn’t terribly informative to humans: Because we’re multiplying (i.e., taking the “product” of) a quantity that depends on X by a quantity that depends on Y and then averaging, you can think of the formula for the covariance as an “average cross product” between and . The covariance has the nice property that, if and are entirely unrelated, then the covariance is exactly zero. If the relationship between them is positive (in the sense shown in Figure 12.4 then the covariance is also positive, and if the relationship is negative then the covariance is also negative. In other words, the covariance captures the basic qualitative idea of correlation. Unfortunately, the raw magnitude of the covariance isn’t easy to interpret as it depends on the units in which and are expressed and, worse yet, the actual units that the covariance itself is expressed in are really weird. For instance, if refers to the dani.sleep variable (units: hours) and refers to the dani.grump variable (units: grumps), then the units for their covariance are . And I have no freaking idea what that would even mean. The Pearson correlation coefficient r fixes this interpretation problem by standardising the covariance, in pretty much the exact same way that the z-score standardises a raw score, by dividing by the standard deviation. However, because we have two variables that contribute to the covariance, the standardisation only works if we divide by both standard deviations. In other words, the correlation between and can be written as follows:
— Just like we saw with the variance and the standard deviation, in practice we divide by rather than .
This is an oversimplification, but it’ll do for our purposes.↩︎Also sometimes written as
where is the slope coefficient and is the intercept (constant) coefficient: ↩︎The
symbol is the Greek letter epsilon. It’s traditional to use or to denote a residual: which in turn means that we can write down the complete linear regression model as: ↩︎Or at least, I’m assuming that it doesn’t help most people. But on the off chance that someone reading this is a proper kung fu master of linear algebra (and to be fair, I always have a few of these people in my intro stats class), it will help you to know that the solution to the estimation problem turns out to be
, where is a vector containing the estimated regression coefficients, is the “design matrix” that contains the predictor variables (plus an additional column containing all ones; strictly is a matrix of the regressors, but I haven’t discussed the distinction yet), and is a vector containing the outcome variable. For everyone else, this isn’t exactly helpful and can be downright scary. However, since quite a few things in linear regression can be written in linear algebra terms, you’ll see a bunch of footnotes like this one in this chapter. If you can follow the maths in them, great. If not, ignore it.↩︎The formula for the general case: The equation that I gave in the main text shows you what a multiple regression model looks like when you include two predictors. Not surprisingly then, if you want more than two predictors all you have to do is add more
terms and more coefficients. In other words, if you have predictor variables in the model then the regression equation look like this: ↩︎And by “sometimes” I mean “almost never”. In practice everyone just calls it “R-squared”.↩︎
The adjusted
value introduces a slight change to the calculation, as follows. For a regression model with predictors, fit to a data set containing observations, the adjusted is: ↩︎Formally, our “null model” corresponds to the fairly trivial “regression” model in which we include 0 predictors and only include the intercept term
: If our regression model has predictors, the “alternative model” is described using the usual formula for a multiple regression model: How can we test these two hypotheses against each other? The trick is to understand that it’s possible to divide up the total variance into the sum of the residual variance and the regression model variance . I’ll skip over the technicalities, since we’ll get to that later when we look at ANOVA in Chapter 13. But just note that And we can convert the sums of squares into mean squares by dividing by the degrees of freedom: So, how many degrees of freedom do we have? As you might expect the df associated with the model is closely tied to the number of predictors that we’ve included. In fact, it turns out that . For the residuals the total degrees of freedom is . Now that we’ve got our mean square values we can calculate an -statistic like this: and the degrees of freedom associated with this are and .↩︎For advanced readers only. The vector of residuals is
. For K predictors plus the intercept, the estimated residual variance is . The estimated covariance matrix of the coefficients is , the main diagonal of which is , our estimated standard errors.↩︎Note that, although jamovi has done multiple tests here, it hasn’t done a Bonferroni correction or anything (see Chapter 13). These are standard one-sample
-tests with a two-sided alternative. If you want to make corrections for multiple tests, you need to do that yourself.↩︎Fortunately, confidence intervals for the regression weights can be constructed in the usual fashion
where is the standard error of the regression coefficient, and is the relevant critical value of the appropriate -distribution. For instance, if it’s a 95% confidence interval that we want, then the critical value is the th quantile of a -distribution with degrees of freedom. In other words, this is basically the same approach to calculating confidence intervals that we’ve used throughout.↩︎Strictly, you standardise all the regressors. That is, every “thing” that has a regression coefficient associated with it in the model. For the regression models that I’ve talked about so far, each predictor variable maps onto exactly one regressor, and vice versa. However, that’s not actually true in general and we’ll see some examples of this later in Chapter 14. But, for now we don’t need to care too much about this distinction.↩︎
Leaving aside the interpretation issues, let’s look at how it’s calculated. What you could do is standardise all the variables yourself and then run a regression, but there’s a much simpler way to do it. As it turns out, the
coefficient for a predictor and outcome has a very simple formula, namely where is the standard deviation of the predictor, and is the standard deviation of the outcome variable . This makes matters a lot simpler.↩︎The way we calculate these is to divide the ordinary residual by an estimate of the (population) standard deviation of these residuals. For technical reasons, the formula for this is:
where in this context is the estimated population standard deviation of the ordinary residuals, and is the “hat value” of the th observation. I haven’t explained hat values to you yet, so this won’t make a lot of sense. For now, it’s enough to interpret the standardised residuals as if we’d converted the ordinary residuals to z-scores.↩︎The formula for doing the calculations this time is subtly different
Notice that our estimate of the standard deviation here is written . What this corresponds to is the estimate of the residual standard deviation that you would have obtained if you just deleted the th observation from the data set. This sounds like the sort of thing that would be a nightmare to calculate, since it seems to be saying that you have to run new regression models (even a modern computer might grumble a bit at that, especially if you’ve got a large data set). Fortunately, this standard deviation estimate is actually given by the following equation: .↩︎ for all values of except . When we just take the natural logarithm (i.e., ln( ).↩︎In jamovi, you can compute this new variable using the formula ‘SQRT(ABS(Residuals))’.↩︎
It’s a bit beyond the scope of this chapter to talk about how to deal with violations of homogeneity of variance, but I’ll give you a quick sense of what you need to consider. The main thing to worry about, if homogeneity of variance is violated, is that the standard error estimates associated with the regression coefficients are no longer entirely reliable, and so your
-tests for the coefficients aren’t quite right either. A simple fix to the problem is to make use of a “heteroscedasticity corrected covariance matrix” when estimating the standard errors. These are often called sandwich estimators, and these can be estimated in (but not directly in jamovi).↩︎The formula for the
-th VIF is: where refers to R-squared value you would get if you ran a regression using as the outcome variable, and all the other variables as the predictors. The idea here is that is a very good measure of the extent to which is correlated with all the other variables in the model. Better yet, the square root of the VIF is pretty interpretable: it tells you how much wider the confidence interval for the corresponding coefficient is, relative to what you would have expected if the predictors are all nice and uncorrelated with one another.↩︎Again, for the linear algebra fanatics: the “hat matrix” is defined to be that matrix
that converts the vector of observed values into a vector of predicted values , such that . The name comes from the fact that this is the matrix that “puts a hat on y”. The hat value of the i-th observation is the i-th diagonal element of this matrix (so technically I should be writing it as rather than ). And here’s how it’s calculated: .↩︎ Notice that this is a multiplication of something that measures the outlier-ness of the observation (the bit on the left), and something that measures the leverage of the observation (the bit on the right).↩︎In jamovi you can save the Cook’s distance values to the data set, then draw a boxplot of the Cook’s distance values to identify the specific outliers. Or you could use a more powerful regression program such as the “car” package in R which has more options for advanced regression diagnostic analysis.↩︎
In the context of a linear regression model (and ignoring terms that don’t depend on the model in any way!), the AIC for a model that has
predictor variables plus an intercept is .↩︎While I’m on this topic I should point out that the empirical evidence suggests that BIC is a better criterion than AIC. In most simulation studies that I’ve seen, BIC does a much better job of selecting the correct model.↩︎
We can fit both models to the data and obtain a residual sum of squares for both models. I’ll denote these as:
and respectively. The superscripting here just indicates which model we’re talking about. Then our statistic is: where is the number of observations, is the number of predictors in the full model (not including the intercept), and is the difference in the number of parameters between the two models. The degrees of freedom here are and . Note that it’s often more convenient to think about the difference between those two values as a sum of squares in its own right. That is: The reason why this is helpful is that we can express as a measure of the extent to which the two models make different predictions about the the outcome variable. Specifically: where is the predicted value for according to model and is the predicted value for according to model .
—
It’s worth noting in passing that this same statistic can be used to test a much broader range of hypotheses than those that I’m mentioning here. Very briefly, notice that the nested model corresponds to the full model when we constrain some of the regression coefficients to zero. It is sometimes useful to construct sub-models by placing other kinds of constraints on the regression coefficients. For instance, maybe two different coefficients might have to sum to zero. You can construct hypothesis tests for those kind of constraints too, but it is somewhat more complicated and the sampling distribution for can end up being something known as the non-central -distribution, which is way beyond the scope of this book! All I want to do is alert you to this possibility.↩︎