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

1. Download week_05_workshop.R from NOW.
2. Move week_05_workshop.R into your RStudio Project folder.
3. Check that the workshop_data folder is also inside your RStudio Project folder.
4. Open your .Rproj file for this module.
5. Open week_05_workshop.R from inside RStudio.

Please wait for each exercise checkpoint before starting that part of the script. Listen to the explanation and instructions first: they help 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.

• \(Y\) is the observed outcome: here, reaction time.
• \(X\) is the observed predictor: here, age.
• Greek letters represent model parameters that we estimate from the data.

• 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?

A simple regression estimates three parameters:

• \(\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

• The intercept is the predicted reaction time when age is 0.
• The age slope means predicted reaction time increases by about 5.78 ms for each additional year of age.

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.

Open week_05_workshop.R and complete Part A1 with me.

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:

The hats mark estimates: \(\widehat{Y}\) is the predicted mean; \(Y\) is the observed outcome.

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 with me.

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 with me.

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.

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.