Linear Regression

Some background information

For this practical, we’ll be using open data from the following source: Wright, D. and London, K. (2011). Modern Regression Techniques Using R, London: Sage. (Chapter 8). This is available here.

This dataset contains information about a randomised controlled trial, in which primary school pupils were asked about the time they spent doing exercises (weekly). Then, their classes (within schools) were divided into 4 different groups and pupils within those classes were given different treatments:

1. Control (no treatment at all)
2. Leaflet (information about exercising)
3. Leaflet+plan (the leaflet and a plan of exercise)
4. Leaflet+plan+quiz (leaflet, plan and a quiz to do about exercising)

Pupils were asked again about the time they spent doing exercises (weekly) some time after the treatments were given.

The variable sqw2 is the outcome (time spent doing exercises after treatment). The variable sqw1 is the measure before the treatment. The variable wcond is the group or condition to which the class of the pupils was assigned

Typical workflow setup

a) Define a working directory

You can use any directory in your computer. As in the example below:

setwd("C:/myfolder")

b) Read in data

This is from Wright and London’s book website companion.

mydata<-read.delim("https://us.sagepub.com/sites/default/files/upm-binaries/26934_exercise.dat")

Note: You can safely ignore the warning when reading in the data.

c) Load packages

You always need to load the packages that you will be using. In this practical, ggplot2 is optional, but we’ll load it anyway:

library(ggplot2)

---

Task 1: Correlation

Before we start any data analysis process, we need to explore our variables. In the case of linear regression, it is good to know how correlated our variables are.

What is the correlation between exercise before and after the intervention?

cor(mydata$sqw2, mydata$sqw1)

## [1] 0.7648664

Correlations go from -1 to 1, so we can safely say that these two variables are highly correlated. This is not surprising at all, since what we’re measuring is behaviour (time spent doing exercises) and you know what they say “the best predictor of future behaviour is past behaviour” ¹.

In any case, this is overly simplistic, so we’ll proceed to run a series of linear regression models.

---

Task 2: Simple linear regression

Simple linear regression has the following algebraic form:

\[Y_i = \beta_0 + \beta_1x_{i1} + \epsilon_i\]

Where:

\(Y_i\) is the dependent variable (outcome of interest)

\(\beta_0\) is the intercept (overall mean)

\(\beta_1\) is the change in the dependent variable when the independent variable changes in one unit

\(x_{i1}\) is the independent variable (aka explanatory variable, predictor or covariate)

\(\epsilon_i\) is the residual

To run a simple linear regression, we need to specify a name for our model (model1) and then call the function lm (linear model). Then we specify the data we will use mydata. Finally, we specify the actual statistical model: sqw2 on the left-hand side as the dependent variable followed by ~ and then the right-hand side contains the independent variable sqw1.

model1 <- lm(data = mydata, sqw2 ~ sqw1)

summary(model1)

## 
## Call:
## lm(formula = sqw2 ~ sqw1, data = mydata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.21789 -0.23311  0.01738  0.19632  1.25155 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.58851    0.04799   12.26   <2e-16 ***
## sqw1         0.71438    0.02688   26.58   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3632 on 501 degrees of freedom
## Multiple R-squared:  0.585,  Adjusted R-squared:  0.5842 
## F-statistic: 706.3 on 1 and 501 DF,  p-value: < 2.2e-16

Question time!

Q1. What is the effect of sqw1 on sqw2?

Q2. How much variance is explained?

---

Task 3: Multiple linear regression

Normally in social sciences, you wouldn’t expect phenomena to be explained by only one variable. In the example we’re working on

\[Y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \dots + \beta_nx_{in} + \epsilon_i\] All the terms mean the same as before, so we’ll focus on the new information:

\(\beta_2\) is the expected increase in the dependent variable conditional on the rest of the variables remaining constant

\(x_{i2}\) is a second independent variable

\(\beta_nx_{in}\) is used to indicate that there can be multiple independent variables

Now, run a multiple linear regression model by adding the experimental condition variable.

Tip: The experimental condition is a categorical variable. In R, categorical variables are referred to as “factor”. To add a factor variable to a regression model, you need to specify it as such:

factor(name_of_variable)

model2 <- lm(data = mydata, sqw2 ~ sqw1 + factor(wcond))

summary(model2)

## 
## Call:
## lm(formula = sqw2 ~ sqw1 + factor(wcond), data = mydata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.29138 -0.20989  0.05553  0.17752  1.19031 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.44117    0.05419   8.142 3.15e-15 ***
## sqw1            0.71615    0.02630  27.225  < 2e-16 ***
## factor(wcond)2  0.15754    0.04431   3.556 0.000413 ***
## factor(wcond)3  0.20734    0.04403   4.709 3.22e-06 ***
## factor(wcond)4  0.21753    0.04525   4.808 2.03e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3536 on 498 degrees of freedom
## Multiple R-squared:  0.6091, Adjusted R-squared:  0.606 
## F-statistic:   194 on 4 and 498 DF,  p-value: < 2.2e-16

Question time!

Q3. What is the effect of condition 2 in relation to the control group?

Q4. What is the effect of condition 3 in relation to the control group?

Q5. What is the effect of condition 4 in relation to the control group?

---

Task 4: Adding interaction effects

Adding interaction is fairly straightforward in R. All you need to do is “multiply” the terms you want to interact (using *).

Take the previous model and add an interaction term between the experimental condition and exercise before the intervention.

model3 <- lm(data = mydata, sqw2 ~ sqw1 + factor(wcond) + sqw1*factor(wcond))

summary(model3)

## 
## Call:
## lm(formula = sqw2 ~ sqw1 + factor(wcond) + sqw1 * factor(wcond), 
##     data = mydata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.29475 -0.20625  0.05175  0.15790  1.27243 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.38351    0.09700   3.954 8.82e-05 ***
## sqw1                 0.75041    0.05460  13.744  < 2e-16 ***
## factor(wcond)2       0.39801    0.13079   3.043  0.00247 ** 
## factor(wcond)3       0.11658    0.13936   0.837  0.40326    
## factor(wcond)4       0.23328    0.13393   1.742  0.08215 .  
## sqw1:factor(wcond)2 -0.13978    0.07225  -1.935  0.05361 .  
## sqw1:factor(wcond)3  0.05949    0.08123   0.732  0.46428    
## sqw1:factor(wcond)4 -0.01006    0.07398  -0.136  0.89184    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3518 on 495 degrees of freedom
## Multiple R-squared:  0.6153, Adjusted R-squared:  0.6099 
## F-statistic: 113.1 on 7 and 495 DF,  p-value: < 2.2e-16

Question time!

Q6. What is the R-squared of this model?

Q7. Does the addition of interactions improve the model?

---

Task 5: Predictions

You can calculate predicted values by using the values of the coefficients. So, for example, you can predict the score of a pupil who scored 2 in the baseline measure and was in the control group as such:

Since the control group is coded 0, the equation for these pupils reduces to:

\({Y}_i = \beta_0 + \beta_1x_{i1}\)

Then

\(\hat{y} = 0.38351 + 0.75041*baseline\)

\(\hat{y} = 0.38351 + 0.75041*2\)

\(\hat{y} = 1.88433\)

The predicted score after the trial for a child in the control group who scored 2 at baseline is then 1.88433.

Question time!

Q8. What is the predicted score for a child in condition 2 that scored 1.5 at baseline?

Plotting the estimated model

To do this, you can retrieve the predictions from your model object and plot them against the observed values of x

mydata$predicted <- fitted(model3)

Remember: You need to load ggplot2 for this!

library(ggplot2)

Then draw a line plot. The resulting line is the predicted effect of sqw1 on sqw2, conditional upon the experimental condition and the interaction between sqw1 and the experimental condition.

plot1 <- ggplot(mydata, aes(y = predicted, x = sqw1)) + 
  geom_smooth(method = "lm", aes(y = predicted, x = sqw1), se=F)
plot1

---

Task 6: Checking assumptions

We are going to focus on the two most basic assumption checks that you can do, which are to do with the residuals. We will use graphic methods, since they are quite straightforward and easily understood.

6.1. Normality of residuals

A key assumption in linear regression is that our residuals are normally distributed. This means roughly that whatever is left unexplained by our model could be thought of as random variation or “white noise”.

This is a very useful gif (found via Giphy) that shows random white noise

The normality of the residuals can be checked by running a formal statistical test, e.g. the Kolmogorov-Smirnov test, but those methods can be overly sensitive to small departures. Graphic methods, such as drawing a histogram or a Q-Q plot, are widely used instead.

To draw a Q-Q plot from your model results, you can quickly use the base package as such:

plot(model3,2)

The result is a Q-Q plot or normal quantile plot. It compares the residuals against the normal distribution. Put simply, the residuals (dots on the plot) should all lie on the diagonal line. Any obvious patterns would indicate that our model does not meet the normality assumption.

In the plot above, we can see that the extremes of the distribution depart from normality, which could be indicating that there are outliers. It may also be the case that we need to control for other variables.

6.2. Homoscedasticity

Homoscedasticity means that residuals should distribute evenly in terms of spread across the predicted values of our model. In other words, if we plot residuals vs fitted values, it should look completely random, i.e. a random cloud of dots.

To check this, we can run the following code:

plot(model3,1)

What we are interested to see here is no obvious patterns. In this specific plot, there seems to be some non-random pattern (not a huge one). The left-hand side of the distribution (lower predicted values) tends to have more positive residuals, whereas the right-hand side has more negative residuals. How do you address this? Usually by adding more variables to the model.

Bonus track

You can add the residuals to your data (as a new variable) to be able to work with them alongside the observed data more easily. To retrieve the residuals of your model, you can do the following:

mydata$residuals <- residuals(model3)

Then you can manipulate them as you wish and, for example, draw plots with ggplot2:

This is to draw a histogram of the residuals with the ggplot2 package:

hist <- ggplot(mydata, aes(residuals)) + geom_histogram()

hist

And this is to check for homoscedasticity:

homoscedasticity <- ggplot(mydata, aes(y = residuals, x = predicted)) + geom_point()

homoscedasticity

Acknowledgement

This practical reproduces some of the materials created by Ana Morales-Gomez, UKDS, The University of Manchester. Those materials are available here.

Thank you ever so much, Ana! (gif via Giphy)

---

1. Sorry, I couldn’t find a precise citation for this. There seems to be a lot of attribution uncertainty regarding this quote.↩︎

 

---

Patricio Troncoso

patricio.troncoso@manchester.ac.uk