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.

Get set up

Open RStudio.
Create a new .Rmd file for this week’s exercises.
Save it somewhere you can find it again.
Give it a clear name (for example, dapr2_lab02.Rmd).
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)

More 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”).

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.

Produce two scatterplots, showing how each of the independent (predictor) variables is individually associated with the outcome variable. Use geom_smooth(method = 'lm', se = FALSE) to add the line of best fit to each scatterplot.

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

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

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.

What kind of correlation does this coefficient represent?

Code hint

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

cor(variable1, 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.

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?

🗂️ See Multiple regression > Model specification flash card.

\[ \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.

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.

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:

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.

	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
R² / R² 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.)

Based on this plot and the relevant model coefficient you identified above, how would you respond to the RQ?

Use the model to make predictions

Question 9

Take the coefficients estimated in mdl, rounded to two decimal points, and substitute them appropriately into the mathematical model expression you wrote in Q4.

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

\[ \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 the wellbeing scores of the following three people (remember you can use the R console as a calculator):

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?

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

For Shuting:

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

Write:

\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.