02: Multiple regression

This week, you’ll add another predictor into the same model you fit last week. You’ll practice interpreting coefficients in this multiple regression model (i.e., a model with more than one predictor), and you’ll learn how to make a nice plot of one of the resulting model-estimated associations.

NoteGet set up

1. Open RStudio.
2. Create a new R Markdown (.Rmd) file for this week’s exercises.
3. Save it somewhere you can find it again.
4. Give it a clear name (for example, dapr2_lab02.Rmd).
5. In the first code chunk, load the packages you’ll need this week (and install them if you don’t have them already):
   • tidyverse
   • patchwork
   • sjPlot

---

Research question (RQ): Is there an association between wellbeing and time spent outdoors, when controlling for the effect that social interaction has on wellbeing?

Data dictionary:

variable	description
age	Age in years of respondent
outdoor_time	Self report estimated number of hours per week spent outdoors
social_int	Self report estimated number of social interactions per week (both online and in-person)
routine	Binary 1=Yes/0=No response to the question 'Do you follow a daily routine throughout the week?'
wellbeing	Warwick-Edinburgh Mental Wellbeing Scale (WEMWBS), a self-report measure of mental health and wellbeing. The scale is scored by summing responses to each item, with items answered on a 1 to 5 Likert scale. The minimum scale score is 14 and the maximum is 70
location	Location of primary residence (City, Suburb, Rural)
steps_k	Average weekly number of steps in thousands (as given by activity tracker if available)

NoteMore detail about this dataset

From the Edinburgh & Lothians, 100 city/suburb residences and 100 rural residences were chosen at random and contacted to participate in the study. The Warwick-Edinburgh Mental Wellbeing Scale (WEMWBS) was used to measure mental health and wellbeing.

Participants filled out a questionnaire including items concerning: estimated average number of hours spent outdoors each week, estimated average number of social interactions each week (whether on-line or in-person), whether a daily routine is followed (yes/no). For those respondents who had an activity tracker app or smart watch, they were asked to provide their average weekly number of steps.

Read in and explore data

Question 1

Read in the data from https://uoepsy.github.io/data/wellbeing_rural.csv and store it in a variable named mwdata (“mw” stands for “mental wellbeing”).

🗂️ See Multiple regression > Example flash card.

Solution

Solution 1.

mwdata <- read_csv("https://uoepsy.github.io/data/wellbeing_rural.csv")

Question 2

Recall our RQ: Is there an association between wellbeing and time spent outdoors, when controlling for the effect that social interaction has on wellbeing?

Identify the three relevant variables (that is, the three columns in mwdata) which we’ll use to address this question. Write down which one is the outcome variable and which two are the predictors.

Create two scatterplots, showing how each of the independent (predictor) variables is individually associated with the outcome variable.

Add on geom_smooth(method = 'lm', se = FALSE) to each plot. This will add the straight line of best fit to each scatterplot.

Bonus challenge: use patchwork to create a single graphic that shows the two scatterplots beside each other.

TipCode hint: patchwork

To patchwork two plots together, you’ll first need to store each of the plots under their own name, for example:

# the first plot
plot_petallength <- iris |>
  ggplot(aes(x = Sepal.Width, y = Petal.Length)) +
  geom_point()

# the second plot
plot_petalwidth <- iris |>
  ggplot(aes(x = Sepal.Width, y = Petal.Width)) +
  geom_point()

Then as long as you’ve got the patchwork library loaded, you can write the two names with a + sign in between, and you’ll get the two plots placed side by side:

# the first plot + the second plot
plot_petallength + plot_petalwidth

(If patchwork is not loaded, you’ll get an error about bad use of arithmetic operator, because R won’t know how to handle the + sign.)

🗂️ See Exploring data visually > Bivariate associations flash card.

Solution

Solution 2. The relevant variables in mwdata:

• wellbeing: outcome
• outdoor_time: one predictor
• social_int: another predictor (included so that we can hold it constant)

wellbeing_outdoor <- 
  ggplot(data = mwdata, aes(x = outdoor_time, y = wellbeing)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "Time spent outdoors \nper week (hours)", y = "Wellbeing score (WEMWBS)")

wellbeing_social <- 
  ggplot(data = mwdata, aes(x = social_int, y = wellbeing)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "Number of social interactions \nper week", y = "Wellbeing score (WEMWBS)")

# place plots adjacent to one another using notation from the `patchwork` package
wellbeing_outdoor + wellbeing_social

Note: You’ll probably have noticed that geom_smooth() makes R print out some extra plotting info, namely:

geom_smooth() using formula = ‘y ~ x’

If you want to suppress that message, write message = FALSE in the header of your code chunk.

Question 3

When we include multiple predictors in a linear model, we want to make sure that they’re aren’t too highly correlated. (In Week 8, we’ll go into more detail about why.)

Use the cor() function to find the correlation coefficient for the correlation between the two predictor variables.

Write down what kind of correlation this coefficient represents. (Positive or negative? Weak, medium, or strong?) If in doubt, check the flash cards.

TipCode hint: cor()

Here’s the basic usage of cor() for getting the correlation between two variables:

cor(variable1, variable2)

If you want to pull each of the two variables out of a dataframe called df:

cor(df$variable1, df$variable2)

Or, if the two variables are the only two columns in a dataframe called df:

cor(df)

🗂️ See Exploring data numerically > Correlation flash card.

Solution

Solution 3. There are a couple ways to get the correlation coefficient for the correlation between outdoor_time and social_int.

You might use the dollar sign notation to extract each column from mwdata:

cor(mwdata$outdoor_time, mwdata$social_int)

[1] -0.03610859

You might also subset mwdata using select() and then pipe the resulting two columns into cor():

mwdata |>
  select(outdoor_time, social_int) |>
  cor()

             outdoor_time  social_int
outdoor_time   1.00000000 -0.03610859
social_int    -0.03610859  1.00000000

(This option gives you a correlation matrix. The 1s on the diagonal indicate that each variable is perfectly correlated with itself, and the off-diagonals represent the correlations between different variables.)

The correlation coefficient is –0.04, which is a small/weak negative correlation: as outdoor time increases, social interaction tends to decrease (or vice versa), but only a tiny bit.

Set up and fit linear model

Question 4

Write the mathematical model formulation for the linear model that will address our RQ.

Which \(\beta\) coefficient will we interpret in order to address our RQ about the relationship between outdoor time and wellbeing, when controlling for social interaction? Write this down.

🗂️ See Multiple regression > Model specification flash card.

Solution

Solution 4. \[ \text{wellbeing} = \beta_0 + (\beta_1 \cdot \text{outdoor\_time}) + (\beta_2 \cdot \text{social\_int}) + \epsilon \]

Write:

$$
\text{wellbeing} = \beta_0 + 
(\beta_1 \cdot \text{outdoor\_time}) +
(\beta_2 \cdot \text{social\_int}) +
\epsilon
$$

The \(\beta\) coefficient that we’ll use to address our RQ is \(\beta_1\), which is the slope of the line that associates outdoor time with wellbeing while holding social interaction constant.

Question 5

Use the function lm() to fit this linear model. Name the result mdl.

🗂️ See Multiple regression > Fitting model flash card.

Solution

Solution 5.

mdl <- lm(wellbeing ~ outdoor_time + social_int, data = mwdata)

Look at model estimates

Question 6

Look at the model summary using summary().

Write one sentence per \(\beta\) coefficient, interpreting what each coefficient means in the context of this data and RQ.

🗂️ See Multiple regression > Interpreting results flash card.

Solution

Solution 6.

summary(mdl)

Call:
lm(formula = wellbeing ~ outdoor_time + social_int, data = mwdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.7611  -3.1308  -0.4213   3.3126  18.8406 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  28.62018    1.48786  19.236  < 2e-16 ***
outdoor_time  0.19909    0.05060   3.935 0.000115 ***
social_int    0.33488    0.08929   3.751 0.000232 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.065 on 197 degrees of freedom
Multiple R-squared:  0.1265,    Adjusted R-squared:  0.1176 
F-statistic: 14.26 on 2 and 197 DF,  p-value: 1.644e-06

• Intercept (\(\beta_0\)): The estimated wellbeing score of somebody with zero hours of outdoor time and zero social interactions per week is 28.62 points.
• Slope over outdoor_time (\(\beta_1\)): Increasing one hour of outdoor time per week is associated with an increase in wellbeing score of 0.20 points, holding the number of weekly social interactions constant.
• Slope over social_int (\(\beta_2\)): Increasing one social interaction per week is associated with an increase in wellbeing score of 0.33 points, holding the number of outdoor hours constant.

(Remember to name the units after you’ve given the number! That is, write “28.62 points” and not just “28.62”.)

Question 7

Now you’ll use a function called tab_model() to create a pretty-printed table of coefficient estimates. It comes from the package sjPlot.

Run the following code and observe what it gives you:

tab_model(mdl)

Caution: tab_model() is very nice, but if you include it in an Rmd file which you then want to knit to PDF, it will cause knitr to crash. So it’s best to use tab_model() only in Rmd files that you want to knit to HTML. If you need to knit to PDF, then you can recreate the output of tab_model() manually.

Solution

Solution 7.

tab_model(mdl)

	wellbeing
Predictors	Estimates	CI	p
(Intercept)	28.62	25.69 – 31.55	<0.001
outdoor time	0.20	0.10 – 0.30	<0.001
social int	0.33	0.16 – 0.51	<0.001
Observations	200
R2 / R2 adjusted	0.126 / 0.118

Question 8

Here you’ll try out a new function called plot_model(), also from the package sjPlot.

Run the following code to generate a plot showing the association between outdoor_time and wellbeing, for a few different values of social_int.

plot_model(
    mdl, 
    type = 'eff', 
    terms = c(
        'outdoor_time', 
        'social_int'
    ),
    show.data = TRUE
)

(If you’re not sure what one of the lines of this code does, then comment it out and run the rest of the code again. Pay attention to what changes.)

Write down brief responses to each of the following questions:

1. Based on the plot, is the effect of outdoor_time the same for all values of social_int? (Hint: A variable’s “effect” is shown by its slope.)
2. Based on the coefficient estimates and/or the plot, is the RQ supported by this analysis?

Solution

Solution 8.

plot_model(
    mdl, 
    type = 'eff', 
    terms = c(
        'outdoor_time', 
        'social_int'
    ),
    show.data = TRUE
)

1. Yes, the effect of outdoor_time is the same for all values of social_int. We can tell because all three coloured lines, which show different values of social_int (the mean in blue, the mean + 1 SD in red, and the mean – 1 SD in green), have the same slope over outdoor_time.

2. This plot shows a positive association between outdoor time and wellbeing, even when controlling for social interaction: no matter how much social interaction people get, people with more outdoor time are estimated to have higher wellbeing scores. And the \(\beta\) coefficient that quantifies the slope of this line is also positive at 0.20. Therefore I’d conclude that yes, there seems to be a positive association between wellbeing and time spent outdoors, controlling for the effect that social interaction also has on wellbeing. (Next week, we’ll learn how to report the corresponding significance test!)

Use the model to make predictions

Question 9

Take the coefficients estimated in mdl and round them to two decimal places.

Substitute them appropriately into the mathematical model expression you wrote in Q4.

🗂️ See Compute model-predicted values > Example flash card.

Solution

Solution 9. To let R do the rounding for you:

coef(mdl) |> round(2)

 (Intercept) outdoor_time   social_int 
       28.62         0.20         0.33 

\[ \text{wellbeing} = 28.62 + (0.20 \cdot \text{outdoor\_time}) + (0.33 \cdot \text{social\_int}) + \epsilon \]

Write:

$$
\text{wellbeing} = 28.62 + 
(0.20 \cdot \text{outdoor\_time}) +
(0.33 \cdot \text{social\_int}) +
\epsilon
$$

Question 10

Last week, you used your linear expression to reconstruct the estimated wellbeing score of somebody whose data we observed (the person whose data was in the third row of mwdata).

Finally, this week, you’ll use the linear expression you wrote in Q9 to predict wellbeing scores for people that we haven’t observed any data from.

Calculate by hand (not using predict()) the wellbeing scores of the following three people:

• Shuting:
  • Outdoor time = 36 hours
  • Social interactions = 19
• Fatima:
  • Outdoor time = 20 hours
  • Social interactions = 15
• Donna:
  • Outdoor time = 1 hour
  • Social interactions = 7

Who has the highest predicted wellbeing score? Who has the lowest? Write your responses down.

🗂️ See Compute model-predicted values > Example flash card.

Solution

Solution 10. For Shuting:

\[ \begin{align} \text{wellbeing}_{\text{Shuting}} &= 28.62 + (0.20 \cdot 36) + (0.33 \cdot 19) \\ &= 42.09 \\ \end{align} \]

If you wanted to present the working-out nicely (which is not required), here’s how you might write it:

\begin{align}
\text{wellbeing}_{\text{Shuting}} &= 28.62 + (0.20 \cdot 36) + (0.33 \cdot 19) \\
  &= 42.09 \\
\end{align}

For Fatima:

\[ \begin{align} \text{wellbeing}_{\text{Fatima}} &= 28.62 + (0.20 \cdot 20) + (0.33 \cdot 15) \\ &= 37.57 \\ \end{align} \]

For Donna:

\[ \begin{align} \text{wellbeing}_{\text{Donna}} &= 28.62 + (0.20 \cdot 1) + (0.33 \cdot 7) \\ &= 31.13 \\ \end{align} \]

Shuting has the highest predicted wellbeing score at 42.09, and Donna has the lowest at 31.13.

Optional: Knit your report to PDF

Note: During the labs, please spend your time on the exercises above, not on rendering your document. This guidance is here in case you have extra time at the end.

Rmarkdown is a useful format because it can be used as a foundation for generating different kinds of documents. For the group report, you’ll be asked to render your Rmd file to PDF. Why not practice now?

First, make sure you have tinytex installed in R. This is what lets you convert things into PDF format.

Run the following code to install tinytex (you’ll only need to run it once).

install.packages("tinytex")
tinytex::install_tinytex()

Then ensure that the “yaml” (rhymes with “mammal”) header—the bit at the top of your Rmd document—looks something like this:

---
title: "this is my report title"
author: "B1234506"
date: "07/09/2025"
output: bookdown::pdf_document2
---

Then follow the instructions from the Rmd bootcamp to knit your Rmd document to PDF.

If you encounter errors, refer to the DAPR1 formatting resources / knitting checklist.

TipWhat to do if you cannot knit to PDF

If you are having issues knitting directly to PDF, try the following:

• Knit to HTML file
  
• Open your HTML in a web-browser (e.g. Chrome, Firefox)
  
• Print to PDF (Ctrl+P, then choose to save to PDF)
  
• Open file to check formatting

TipHiding code and/or output from code chunks

Review Lesson 5 of the rmd bootcamp for a detailed description/worked examples.

Hiding R Code

To not show the code of an R code chunk, and only show the output, write:

```{r, echo=FALSE}
# code goes here
```

Hiding R Output:

To show the code of an R code chunk, but hide the output, write:

```{r, results='hide'}
# code goes here
```

Hiding R code AND output

To hide both code and output of an R code chunk, write:

```{r, include=FALSE}
# code goes here
```