1. The Intuition

Regression asks: if I know x, what is my best guess for y? The word simple means one predictor; linear means the relationship is a straight line.

Environmental applications are everywhere. Does PM2.5 increase as traffic volume rises? Does soil carbon decrease with distance from a riparian corridor? Does tree height predict biomass? In every case the method is the same: fit a line that minimizes prediction error, then ask how confident we are in its slope.

# Ozone vs. solar radiation from the built-in airquality dataset
data(airquality)
aq <- na.omit(airquality[, c("Ozone", "Solar.R", "Temp", "Wind")])
head(aq)

##   Ozone Solar.R Temp Wind
## 1    41     190   67  7.4
## 2    36     118   72  8.0
## 3    12     149   74 12.6
## 4    18     313   62 11.5
## 7    23     299   65  8.6
## 8    19      99   59 13.8

Ozone concentration vs. temperature in New York, summer 1973. Each point is one day. The dashed line is a visual guess; the solid line is the least-squares fit.

2. The Model

The simple linear regression model is:

\[y_i = \beta_0 + \beta_1 x_i + \varepsilon_i\]

where:

• \(y_i\) is the observed response for observation \(i\)
• \(x_i\) is the predictor value
• \(\beta_0\) is the intercept — the predicted value of y when x = 0
• \(\beta_1\) is the slope — the change in y for a one-unit increase in x
• \(\varepsilon_i \sim N(0, \sigma^2)\) is the error term — the part of y not explained by x

The ordinary least-squares (OLS) estimator chooses \(\hat\beta_0\) and \(\hat\beta_1\) to minimize the sum of squared residuals:

\[\text{SSR} = \sum_{i=1}^n (y_i - \hat y_i)^2 \quad \text{where} \quad \hat y_i = \hat\beta_0 + \hat\beta_1 x_i\]

The closed-form solution is:

\[\hat\beta_1 = \frac{\sum(x_i - \bar x)(y_i - \bar y)}{\sum(x_i - \bar x)^2} = \frac{S_{xy}}{S_{xx}}\]

\[\hat\beta_0 = \bar y - \hat\beta_1 \bar x\]

x <- aq$Temp
y <- aq$Ozone

b1 <- sum((x - mean(x)) * (y - mean(y))) / sum((x - mean(x))^2)
b0 <- mean(y) - b1 * mean(x)

knitr::kable(
  data.frame(
    Coefficient = c("Slope (β₁)", "Intercept (β₀)"),
    Estimate    = round(c(b1, b0), 3),
    Interpretation = c("ppb per °F", "ppb when Temp = 0°F")
  ),
  caption = "OLS estimates computed by hand."
)

3. Fitting in R with lm()

In practice, always use lm(). The formula syntax y ~ x maps directly onto the statistical model \(y = \beta_0 + \beta_1 x + \varepsilon\).

mod <- lm(Ozone ~ Temp, data = aq)

The fitted object mod contains everything: coefficients, residuals, fitted values, the variance–covariance matrix of the estimates, and more.

coef(mod)          # intercept and slope

## (Intercept)        Temp 
##  -147.64607     2.43911

fitted(mod)[1:6]   # predicted values yhat

##         1         2         3         4         7         8 
## 15.774291 27.969841 32.848061  3.578742 10.896071 -3.738588

residuals(mod)[1:6] # observed - predicted

##          1          2          3          4          7          8 
##  25.225709   8.030159 -20.848061  14.421258  12.103929  22.738588

4. Reading summary() Output

summary(mod)

## 
## Call:
## lm(formula = Ozone ~ Temp, data = aq)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.922 -17.459  -0.874  10.444 118.078 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -147.6461    18.7553  -7.872 2.76e-12 ***
## Temp           2.4391     0.2393  10.192  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.92 on 109 degrees of freedom
## Multiple R-squared:  0.488,  Adjusted R-squared:  0.4833 
## F-statistic: 103.9 on 1 and 109 DF,  p-value: < 2.2e-16

The output has four important sections:

Residuals. A five-number summary of the residuals. Ideally symmetric around zero. Large asymmetry suggests the normality assumption is violated.

Coefficients table. Each row is one model parameter.

Residual standard error (RSE). \(\hat\sigma\) — typical size of a residual. Smaller is better; units match y.

R-squared. The proportion of variance in y explained by x. \(R^2 = 1 - \text{SS}_\text{res} / \text{SS}_\text{tot}\). Range 0–1; higher is better, but high \(R^2\) does not guarantee the model is correct.

F-statistic. Tests whether the overall model explains more variance than a model with no predictors (intercept only). For simple regression, \(F = t^2\).

5. The Four OLS Assumptions

Linear regression produces valid inference only when four assumptions hold. They are easiest to remember with the acronym LINE:

The most important are linearity and equal variance; the t-tests and confidence intervals are fairly robust to modest normality violations when n is large.

6. Residual Diagnostics

plot(mod) produces four standard diagnostic plots. Always examine all four.

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(mod)

Four standard regression diagnostic plots for the ozone model. (1) Residuals vs. Fitted checks linearity and equal variance. (2) Q-Q plot checks normality. (3) Scale-Location is a second variance check. (4) Residuals vs. Leverage identifies influential points.

par(mfrow = c(1, 1))

# Cook's distance > 4/n is a common threshold
n  <- nrow(aq)
cd <- cooks.distance(mod)
influential <- which(cd > 4 / n)

if (length(influential) > 0) {
  knitr::kable(
    aq[influential, ],
    caption = paste0("Observations with Cook's distance > 4/", n,
                     " (potentially influential).")
  )
}

7. Inference: Confidence Intervals and the F-test

confint(mod, level = 0.95)

##                   2.5 %      97.5 %
## (Intercept) -184.818372 -110.473773
## Temp           1.964787    2.913433

Visualizing the fitted line with uncertainty band

new_x  <- data.frame(Temp = seq(min(aq$Temp), max(aq$Temp), length.out = 200))
ci_fit <- predict(mod, newdata = new_x, interval = "confidence")
pi_fit <- predict(mod, newdata = new_x, interval = "prediction")

par(las = 1, mar = c(4.5, 4.5, 2, 1))
plot(Ozone ~ Temp, data = aq, pch = 20,
     col = rgb(0.29, 0.49, 0.37, 0.55),
     xlab = "Temperature (°F)", ylab = "Ozone (ppb)",
     main = "Fitted regression with CI and PI")

# Prediction interval ribbon
polygon(c(new_x$Temp, rev(new_x$Temp)),
        c(pi_fit[, "lwr"], rev(pi_fit[, "upr"])),
        col = rgb(0.2, 0.35, 0.25, 0.08), border = NA)

# Confidence interval ribbon
polygon(c(new_x$Temp, rev(new_x$Temp)),
        c(ci_fit[, "lwr"], rev(ci_fit[, "upr"])),
        col = rgb(0.2, 0.35, 0.25, 0.20), border = NA)

lines(new_x$Temp, ci_fit[, "fit"], col = "#1e3a2f", lwd = 2)
lines(new_x$Temp, pi_fit[, "lwr"], col = "#4a7c5e", lwd = 1, lty = 2)
lines(new_x$Temp, pi_fit[, "upr"], col = "#4a7c5e", lwd = 1, lty = 2)

legend("topleft",
       legend = c("Fitted line", "95% CI (mean)", "95% PI (individual)"),
       col = c("#1e3a2f", rgb(0.2,.35,.25,.5), "#4a7c5e"),
       lty = c(1, 1, 2), lwd = c(2, 6, 1), bty = "n", cex = 0.85)

Fitted regression line (dark) with 95% confidence band (grey) and 95% prediction interval (dashed). The CI captures uncertainty in the mean response; the PI captures uncertainty in a single future observation.

📌 CI vs. PI — a critical distinction The confidence interval captures uncertainty about the mean response at a given x (the population average). The prediction interval captures uncertainty about a single new observation. The PI is always wider because it must account for both the uncertainty in the mean and the random scatter of individual observations around that mean.

---

8. Making Predictions

# Predict ozone at specific temperatures
predict_data <- data.frame(Temp = c(65, 75, 85, 95))
preds <- predict(mod, newdata = predict_data, interval = "confidence")

result <- cbind(predict_data, round(preds, 1))
names(result) <- c("Temp (°F)", "Predicted Ozone (ppb)", "95% CI Lower", "95% CI Upper")
knitr::kable(result, caption = "Predicted mean ozone at selected temperatures.")

# 1. Explore
plot(y ~ x, data = mydata)

# 2. Fit
mod <- lm(y ~ x, data = mydata)

# 3. Summarize
summary(mod)
confint(mod)

# 4. Diagnose
par(mfrow = c(2, 2)); plot(mod); par(mfrow = c(1, 1))

# 5. Predict
predict(mod, newdata = data.frame(x = new_values), interval = "confidence")

mod <- lm(Ozone ~ Temp, data = airquality)
summary(mod)
confint(mod)

par(mfrow = c(2,2))
plot(mod)