• Connect the lecture idea of a line of best fit to R code.
• Fit a simple linear regression with lm().
• Extract coefficients with coef() and residual SD with sigma().
• Use coefficients to calculate predicted means by hand.
• Use predict() for simple and multiple regression predictions.

Check the online lecture before the workshop. Complete the Week 5 problem set after this workshop.

The lecture introduces regression as a normal linear model.

• The line of best fit minimises the sum of squared residuals.
• The intercept is the predicted outcome when the predictor is 0.
• The slope is the expected outcome change for a one-unit predictor change.
• The residual standard deviation describes variability around the regression line.
• With multiple predictors, each slope describes one predictor while the other predictors are held constant.

Workshop focus: fit the model, extract the pieces, and use them for prediction.

For each regression model, use the same workflow:

1. identify the outcome and predictor variables.
2. fit the model with lm().
3. inspect coefficients with the printed model and coef().
4. inspect residual SD with sigma().
5. calculate predictions manually and with predict().

The Week 5 problem set asks for exactly these operations.

Work from inside your RStudio Project folder.

Open your .Rproj file for this module and check that workshop_data is inside the project folder.

1. Download week_05_workshop.R from NOW.
2. Move it into your RStudio Project folder using your file explorer.
3. Open the script in RStudio.

Please wait for each exercise checkpoint before starting that part of the script. The explanation helps you understand what the exercise is meant to show and what to look for in the output.

Run one line at a time and check the console after each line.

Part A of week_05_workshop.R follows the slides. We will complete it together with the Blomkvist data.

blomkvist <- read_csv("workshop_data/analysis_blomkvist.csv")
blomkvist <- drop_na(blomkvist, rt, age, medicine)

How does reaction time change with age?

The Blomkvist data are adapted from a reaction-time and aging study using the Nintendo Wii Balance Board (Blomkvist et al., 2017).

For simple linear regression in R, the formula has this structure:

outcome ~ predictor

For the Blomkvist example:

rt ~ age

• rt is the numeric outcome.
• age is the numeric predictor.
• The outcome goes on the left of ~; the predictor goes on the right.

m_age <- lm(rt ~ age, data = blomkvist)
m_age

lm() fits the linear model. Printing m_age shows the estimated regression coefficients.

This week is about the estimated model parameters and how to use them for prediction.

Complete the data-loading steps and Part A1 in week_05_workshop.R with me.

After fitting and printing m_age, extract its coefficients and residual SD:

b_age <- coef(m_age)
b_age
sigma(m_age)

• Intercept: predicted mean when the predictor is 0.
• Slope: change in the predicted mean per one-unit increase in the predictor.
• Sigma: spread of observations around the line.

coef(m_age)
# (Intercept)        age
#  314.260054   5.778301

sigma(m_age)
# 149.2607

Discuss in pairs / small groups first, before sharing the answers with me:

• Where is the intercept? What does the upward slope tell us?
• What does the ribbon tell us about reaction times?

Greek letters represent model parameters that we estimate from the data.

• \(\beta_0\): the intercept.
• \(\beta_1\): the slope for the predictor.
• \(\sigma\): the residual standard deviation.

coef() extracts the estimated betas.

coef(m_age)
# (Intercept)         age
#  314.260054    5.778301

• Intercept: predicted reaction time at age 0.
• Age slope: about 5.78 ms longer predicted reaction time per additional year.

sigma() extracts the estimated residual standard deviation.

sigma(m_age)
# 149.2607

For people of the same age, reaction times vary around the predicted mean.

Here, sigma is about 149 ms: the differences between observed reaction times and the model’s predicted means have an estimated standard deviation of 149 ms.

In this model, the predicted mean changes with age, but the residual SD is the same at each age.

b_age <- coef(m_age)

b_age contains two numbers.

b_age[1]
# 314.2601

b_age[2]
# 5.778301

Read this as:

• b_age[1]: intercept
• b_age[2]: slope for age

The predicted mean is:

• \(Y\) is the observed outcome: here, reaction time.
• \(X\) is the observed predictor: here, age.
• The hats mark estimates: \(\widehat{Y}\) is the predicted mean.

In R, for age 60:

b_age[1] + b_age[2] * 60
# 660.9581

This is the predicted mean reaction time for participants aged 60.

Start at the intercept, then add the slope multiplied by age:

\[ \widehat{Y} = \widehat{\beta}_0 + \widehat{\beta}_1 \times 60 \approx 660.958\text{ ms} \]

We want the predicted mean reaction time at age 60.

1. Create a small data frame using tibble().
2. Name its column age, exactly as in the model formula.
3. Enter the age for which you want a prediction, then inspect the data frame.

new_age <- tibble(age = 60)
new_age
# A tibble: 1 x 1
#     age
#   <dbl>
# 1    60

new_age is a name we choose for the data frame. The column name age must match the model. We do not supply rt: that is what we want to predict.

predict(m_age, newdata = new_age)
# 660.9581

• m_age tells R which fitted model to use.
• newdata = new_age supplies the predictor values in our new data frame.
• R uses the estimated intercept and slope to calculate the predicted mean.

The answer matches the manual calculation: about 660.958 ms at age 60.

Continue Part A2 in week_05_workshop.R.

Take a short break before we continue.

After the break, we will add a second predictor, then practise the same workflow independently.

For more than one predictor, add predictors with +.

rt ~ age + medicine

The fitted equation is:

Each coefficient describes one predictor while the other predictor is in the model.

Here, \(Y\) is reaction time, \(X_1\) is age, and \(X_2\) is medicine. We estimate four parameters: \(\beta_0\), \(\beta_1\), \(\beta_2\), and \(\sigma\).

m_age_meds <- lm(rt ~ age + medicine, data = blomkvist)
m_age_meds

b_age_meds <- coef(m_age_meds)
b_age_meds
# (Intercept)         age    medicine
#  323.128523    5.421639    6.284839

sigma(m_age_meds)
# 148.9805

For age 60 and medicine 3:

b_age_meds[1] + b_age_meds[2] * 60 + b_age_meds[3] * 3
# 667.2814

With two predictors, the new data frame must contain both predictor columns, named exactly as in the model: age and medicine.

new_person <- tibble(age = 60, medicine = 3)
new_person
predict(m_age_meds, newdata = new_person)
# 667.2814

Continue Part A3 in week_05_workshop.R.

Part B asks whether participant accuracy predicts mean reaction time rt.

Recommended reading: Andrews (2021), Chapter 9, Normal Linear Models.

Problem-set reminder: complete the Week 5 problem set after this workshop. Focus on fitting regression models, extracting beta coefficients and sigma, and using the model for predictions.

If you have not yet completed MySay, open Student voice on this module’s NOW page to find the survey. Tell us what has helped your learning, what you have not found helpful, and what you feel would help you learn. Your feedback is anonymous, and every response matters.

Before the problem set, you should be able to:

• fit lm(y ~ x, data = d).
• fit lm(y ~ x1 + x2, data = d).
• extract beta coefficients with coef().
• extract the residual SD parameter with sigma().
• calculate a prediction by hand from coefficients.
• use predict() with a small newdata data frame.