No installation needed

Scan the QR code or open the link below to try Jupyter in your browser.

https://jupyter.org

Learning goals

• Compute and interpret \(R^2\) from a fitted model.
• Compare \(R^2\) and adjusted \(\bar{R}^2\) when adding regressors.
• Fit a multiple linear regression model in R.
• Interpret partial slopes (ceteris paribus).
• See omitted variable bias in action by comparing a short and a long regression.
• Use anova() to inspect the variance decomposition (SST, SSE, SSR).

Pre-loaded R data

We will use airquality, available in base R.

It contains daily air quality measurements in New York, May–September 1973.

data(airquality)

?airquality   # run this in your session to read the documentation

Variables

Our research question:

What predicts ozone concentration? We start with temperature alone, then add more variables.

A note on missing values

airquality has missing values (NA) in Ozone and Solar.R.

lm() handles them automatically by dropping incomplete rows — but it is good practice to check:

# How many complete cases do we have?
sum(complete.cases(airquality))

[1] 111

# Work with complete cases only
aq <- airquality[complete.cases(airquality), ]
nrow(aq)

[1] 111

Explore the data

1. Create a scatter plot of Ozone (Y) against Temp (X).
2. Describe the direction and strength of the relationship.
3. Does a linear model seem reasonable?

ggplot(aq, aes(x = Temp, y = Ozone)) +
  geom_point() +
  theme_minimal()

plot(aq$Temp, aq$Ozone)

Fit a simple linear regression

Fit the model:

\[\text{Ozone}_i = \beta_0 + \beta_1\,\text{Temp}_i + \varepsilon_i\]

Tasks

• Estimate the model with lm().
• Display the full summary().
• Extract the \(R^2\).

slr <- lm(Ozone ~ Temp, data = aq)
summary(slr)
summary(slr)$r.squared

Interpret the output

Using the output from Exercise 2, answer:

• What does the slope \(\hat{\beta}_1\) mean in this context?
• What proportion of the variation in ozone is explained by temperature alone?
• Is the slope statistically significant? How do you know?

coef(slr)
summary(slr)$r.squared
summary(slr)$coefficients   # shows t-statistics and p-values

The variance decomposition

R can report SST, SSE, and SSR directly.

Tasks

1. Use anova(slr) and examine the output.
2. Identify which row corresponds to SSE and which to SSR.
3. Verify manually that \(R^2 = \text{SSE}/\text{SST}\).

anova(slr)

# SSE and SSR are in the "Sum Sq" column
a   <- anova(slr)
SSE <- a["Temp",    "Sum Sq"]
SSR <- a["Residuals", "Sum Sq"]
SST <- SSE + SSR

SSE / SST   # should equal summary(slr)$r.squared

Residual diagnostics

1. Plot residuals against fitted values — does the spread look constant?
2. Create a Q-Q plot of the residuals — do they look approximately normal?
3. What concerns, if any, do you notice?

par(mfrow = c(2, 2))
plot(slr)
par(mfrow = c(1, 1))

plot(fitted(slr), resid(slr))
abline(h = 0, lty = 2)

qqnorm(resid(slr))
qqline(resid(slr))

Add a second regressor

Now fit the multiple regression model including Wind:

\[\text{Ozone}_i = \beta_0 + \beta_1\,\text{Temp}_i + \beta_2\,\text{Wind}_i + \varepsilon_i\]

Tasks

• Fit the model with lm().
• Compare \(R^2\) and \(\bar{R}^2\) with the SLR from Exercise 2.
• How did the coefficient on Temp change? Why?

mlr1 <- lm(Ozone ~ Temp + Wind, data = aq)
summary(mlr1)

# Compare R² and adjusted R²
summary(slr)$r.squared;      summary(slr)$adj.r.squared
summary(mlr1)$r.squared;     summary(mlr1)$adj.r.squared

Add a third regressor

Now include Solar.R as well:

\[\text{Ozone}_i = \beta_0 + \beta_1\,\text{Temp}_i + \beta_2\,\text{Wind}_i + \beta_3\,\text{Solar.R}_i + \varepsilon_i\]

Tasks

• Fit and summarise the model.
• Build a table comparing \(R^2\) and \(\bar{R}^2\) across all three models.
• Which model would you prefer, and why?

mlr2 <- lm(Ozone ~ Temp + Wind + Solar.R, data = aq)
summary(mlr2)

# Comparison table
data.frame(
  Model    = c("SLR (Temp)", "MLR (Temp+Wind)", "MLR (Temp+Wind+Solar.R)"),
  R2       = c(summary(slr)$r.squared,
               summary(mlr1)$r.squared,
               summary(mlr2)$r.squared),
  Adj_R2   = c(summary(slr)$adj.r.squared,
               summary(mlr1)$adj.r.squared,
               summary(mlr2)$adj.r.squared)
)

Omitted variable bias in action

In Exercise 6 you noticed the coefficient on Temp changed when Wind was added.

Tasks

1. Extract \(\hat{\beta}_{\text{Temp}}\) from the SLR (Exercise 2) and from the MLR with Wind (Exercise 6).
2. Are they different? In which direction?
3. Check whether Temp and Wind are correlated — does the sign of the correlation explain the direction of the bias?

coef(slr)["Temp"]
coef(mlr1)["Temp"]

cor(aq$Temp, aq$Wind)

Prediction with MLR

Using your full model from Exercise 7, predict ozone concentration for two new days:

Tasks

• Obtain point predictions.
• Obtain 95% prediction intervals.
• Which day is predicted to have higher ozone? Does it make sense?

new_days <- data.frame(Temp = c(80, 95), Wind = c(10, 5), Solar.R = c(200, 250))
predict(mlr2, new_days, interval = "prediction", level = 0.95)

Core commands

# Load and clean data
data(airquality)
aq <- airquality[complete.cases(airquality), ]

# Simple linear regression
slr  <- lm(Ozone ~ Temp, data = aq)
summary(slr)

# Multiple linear regression
mlr1 <- lm(Ozone ~ Temp + Wind, data = aq)
mlr2 <- lm(Ozone ~ Temp + Wind + Solar.R, data = aq)
summary(mlr2)

# Goodness of fit
summary(slr)$r.squared
summary(slr)$adj.r.squared

# Variance decomposition
anova(slr)

# Coefficients
coef(mlr2)

# Predictions
new_days <- data.frame(Temp = c(80, 95), Wind = c(10, 5), Solar.R = c(200, 250))
predict(mlr2, new_days, interval = "prediction")

Useful checks

# Residual plots
par(mfrow = c(2, 2))
plot(mlr2)
par(mfrow = c(1, 1))

# Manual residual vs fitted
plot(fitted(mlr2), resid(mlr2))
abline(h = 0, lty = 2)

# Normality check
qqnorm(resid(mlr2))
qqline(resid(mlr2))

# Correlation between regressors (multicollinearity check)
cor(aq[, c("Temp", "Wind", "Solar.R")])

# Compare R² across models
data.frame(
  Model  = c("SLR", "MLR (2 vars)", "MLR (3 vars)"),
  R2     = c(summary(slr)$r.squared,
             summary(mlr1)$r.squared,
             summary(mlr2)$r.squared),
  Adj_R2 = c(summary(slr)$adj.r.squared,
             summary(mlr1)$adj.r.squared,
             summary(mlr2)$adj.r.squared)
)

Homework (not graded)

The full airquality dataset also has Month as a variable.

1. Add Month to the MLR model (treat it as numeric for now) and check whether \(\bar{R}^2\) improves.
2. Look at the scatter plot of Ozone vs Temp coloured by Month. Do you notice any pattern?
3. We will discuss in the next lecture whether controlling for Month makes conceptual sense, and what to do when a regressor is categorical rather than continuous.

ggplot(airquality[complete.cases(airquality),], aes(x = Temp, y = Ozone, color = factor(Month))) +
  geom_point(size = 2.5) +
  theme_minimal() +
  labs(color = "Month")