In Part 1, we fitted normal linear models and inspected coefficient tables.
We focused on questions like:
Today we use the same fitted models, but ask a more concrete question:
What does the model predict for meaningful predictor values?
Today we use fitted linear models to make predictions.
By the end, you should be able to:
predict() for meaningful predictor valuesggplot()Coefficients are compact, but they can be abstract.
Predictions show what the model expects for meaningful predictor values.
Good prediction plots help readers see:
From Part 1, a linear model has two connected pieces:
\[ y_i \sim \mathcal{N}(\mu_i, \sigma) \]
\[ \mu_i = \beta_0 + \beta_1 x_i \]
For predictions, we use the fitted model:
Example model:
model_age <- lm(rt ~ age, data = blomkvist) coef(model_age)
(Intercept) age
314.74 5.77
Formula:
\[ \hat{rt}_i = \hat{\beta}_0 + \hat{\beta}_\text{age} \cdot \text{age}_i \]
To calculate a prediction:
Formula:
\[ \hat{rt}_i = \hat{\beta}_0 + \hat{\beta}_\text{age} \cdot \text{age}_i \]
After running a linear model, we found:
Question:
What is the predicted reaction time for a 40-year-old participant?
Working and answer:
Example model:
model_sex <- lm(rt ~ sex, data = blomkvist) coef(model_sex)
(Intercept) sexmale
664.8 -64.9
The predictor sex has two levels: female and male.
For categorical predictors, R creates 0/1 dummy codes. By default, levels are ordered alphabetically, so female is the reference level and male is the comparison level:
Formula:
\[ \hat{rt}_i = \hat{\beta}_0 + \hat{\beta}_\text{sex} \cdot \text{sexmale}_i \]
Formula:
\[ \hat{rt}_i = \hat{\beta}_0 + \hat{\beta}_\text{sex} \cdot \text{sexmale}_i \]
After running a linear model, we found:
Question:
What is the predicted reaction time for a male participant?
Remember: male is coded as \(\text{sexmale}_i = 1\).
Working and answer:
Example model:
model_age_sex <- lm(rt ~ age + sex, data = blomkvist) coef(model_age_sex)
(Intercept) age sexmale
339.62 5.74 -59.29
Formula:
\[ \hat{rt}_i = \hat{\beta}_0 + \hat{\beta}_\text{age} \cdot \text{age}_i + \hat{\beta}_\text{sex} \cdot \text{sexmale}_i \]
The age term changes with age. The sex term adds the group difference for the comparison group.
Formula:
\[ \hat{rt}_i = \hat{\beta}_0 + \hat{\beta}_\text{age} \cdot \text{age}_i + \hat{\beta}_\text{sex} \cdot \text{sexmale}_i \]
After running a linear model, we found:
Question:
What is the predicted reaction time for a 40-year-old male participant?
Remember: male is coded as \(\text{sexmale}_i = 1\).
Working and answer:
For a simple model, predictions come from the coefficients.
coefs <- coef(model_age_sex) coefs
(Intercept) age sexmale
339.62 5.74 -59.29
# Female is the reference group in this data. coefs[1] + coefs["age"] * 40
(Intercept)
569
predict() is useful because it applies the same logic reliably for many rows and more complex models.
Set up a data frame with hypothetical predictor combinations:
newdata <- data.frame(age = c(20, 40, 60),
sex = c("female", "female", "female"))
newdata
age sex 1 20 female 2 40 female 3 60 female
Predicted average outcomes
predict(model_age_sex, newdata = newdata)
1 2 3 454 569 684
Add Confidence Intervals
predict(model_age_sex, newdata = newdata, interval = "confidence")
fit lwr upr 1 454 415 494 2 569 542 596 3 684 661 707
The confidence interval describes uncertainty in the estimated mean.
prediction_grid <- expand.grid(age = seq(20, 80, by = 5),
sex = c("female", "male"))
head(prediction_grid)
age sex 1 20 female 2 25 female 3 30 female 4 35 female 5 40 female 6 45 female
This creates all combinations of age and sex we want to show.
predicted_values <- predict(model_age_sex,
newdata = prediction_grid,
interval = "confidence")
predictions <- data.frame(prediction_grid, predicted_values)
head(predictions)
age sex fit lwr upr 1 20 female 454 415 494 2 25 female 483 447 519 3 30 female 512 479 544 4 35 female 540 511 570 5 40 female 569 542 596 6 45 female 598 573 623
If you were here yesterday, open part-2-predictions-no-interaction.Rmd from the week-02-exercises folder in your RStudio Project.
If you missed yesterday, download week-02-exercises.zip from the NOW learning room, unzip it, and move the folder into your RStudio Project before opening the Rmd file.
You will:
rt ~ age + sexA ggplot is built in layers.
Start with observed data, then add layers with +.
Key parts:
ggplot() starts the plotaes() maps variables to the plotgeom_point() adds observed datageom_line() adds model predictionsdata = predictions and y = fit adds predicted value from predict()ggplot(blomkvist, aes(x = age, y = rt, colour = sex)) +
geom_point() +
geom_line(data = predictions,
aes(y = fit))
The confidence interval is added with a ribbon.
geom_ribbon() needs three y-values from the prediction data:
fit: predicted meanlwr: lower confidence interval limitupr: upper confidence interval limitThe ribbon layer uses inherit.aes = FALSE because it needs its own mapping: ymin = lwr and ymax = upr.
head(predictions)
age sex fit lwr upr 1 20 female 454 415 494 2 25 female 483 447 519 3 30 female 512 479 544 4 35 female 540 511 570 5 40 female 569 542 596 6 45 female 598 573 623
... + # as before
geom_ribbon(data = predictions,
aes(x = age,
ymin = lwr,
ymax = upr,
fill = sex),
inherit.aes = FALSE) 
A model without an interaction uses +:
lm(rt ~ age + sex, data = blomkvist)
This gives separate age and sex terms, but the age slope is the same for both sex groups.
A model with an interaction uses *:
lm(rt ~ age * sex, data = blomkvist)
age * sex means: include age, include sex, and include age:sex.
Interactions are also often called moderation. Combinations of continuous and categorical predictor are also called varying-intercepts and varying-slopes model:
age + sex)age * sex)In this model, age:sexmale is the interaction term. The value is the difference in age slopes between males and females.
Its estimate is -1.344 ms per year.
coef(model_age_sex_interaction)
(Intercept) age sexmale age:sexmale
309.56 6.27 16.05 -1.34
Without an interaction, the model predicts parallel lines.
The sex coefficient shifts the line up or down, but the age slope is shared.
*I’m not showing raw data, CIs, or any other form of dispersion in this plot. Don’t do this; it’s poor practice. I had to do it to emphasise the difference between predictions with and without interaction.
With an interaction, the model can predict different age slopes for the two sex groups.
In this fitted model, the age-by-sex interaction estimate is -1.344 ms per year.
That means the male age slope is -1.344 ms per year different from the female age slope.
*I’m not showing raw data, CIs, or any other form of dispersion in this plot. Don’t do this; it’s poor practice. I had to do it to emphasise the difference between predictions with and without interaction.
Assume female is the reference group. The interaction asks:
How much different is the age effect for males compared to females?
| Interaction estimate | Meaning |
|---|---|
10 |
male age slope is 10 ms per year more positive than the female age slope |
1 |
male age slope is 1 ms per year more positive than the female age slope |
0 |
male and female age slopes are about the same |
-10 |
male age slope is 10 ms per year less positive than the female age slope |
Positive and negative signs describe the difference between slopes, not just which group is higher.
The reference group matters: the ordinary age coefficient is the age slope for the reference group, and the interaction coefficient tells us how much the comparison group slope differs from that reference slope.
Outcome: sleep duration in hours. Predictors: stress from 0 to 10 and caffeine with levels afternoon and evening.
Reference level: caffeine = afternoon; here caffeine_i = 0 for afternoon and 1 for evening.
\[ \widehat{sleep}_i = \hat{\beta}_0 + \hat{\beta}_\text{stress} \cdot \text{stress}_i + \hat{\beta}_\text{caffeine} \cdot \text{caffeine}_i + \hat{\beta}_\text{stress:caffeine} \cdot \text{stress}_i \cdot \text{caffeine}_i \]
After fitting the model, we find:
-0.20-0.10What does the interaction estimate mean?
A. In the evening-caffeine group, the stress slope is -0.10 hours per stress point.
B. In the afternoon-caffeine group, the stress slope is -0.10 hours per stress point more negative than in the evening-caffeine group.
C. In the evening-caffeine group, the stress slope is -0.10 hours per stress point more negative than in the afternoon-caffeine group.
Then open part-2-predictions-interaction.Rmd.
You will:
rt ~ age + sex and rt ~ age * sexToday we used prediction plots to see what different models imply.
This helps us understand an interaction, but it does not by itself answer whether the interaction is statistically justified.
In Week 3, we will talk about model comparison.
Model comparison asks whether adding a term, such as an interaction, improves model fit enough to justify the more complex model.
If you brought your own dataset, open:
part-2-formative-own-data-guidance.RmdUse it to create predictions for your own formative model and draft a short interpretation paragraph.
If you do not yet have your own dataset, open:
part-2-formative-no-own-data.RmdUse the Blomkvist data to practise the same prediction workflow while you keep looking for a dataset.
Reproduce the prediction workflow with your own dataset.
Try to:
predict() to calculate predicted valuesIf you do not yet have your own dataset, complete part-2-formative-no-own-data.Rmd.
Read after the workshop:
Use these to revise why predictions are often clearer than coefficient tables alone.
Andrews, M. (2021). Doing data science in R: An introduction for Social Scientists. SAGE Publications Ltd.
Field, A., Miles, J., & Field, Z. (2012). Discovering statistics using R. Sage publications.
Gelman, A., Hill, J., & Vehtari, A. (2020). Regression and other stories. Cambridge University Press.