🗓️ Week 04
Linear Regression in R: Theory and Practice
STAT 204 – Introduction to Statistical Data Analysis
19 Oct 2025
Session 1: Simple & Multiple Linear Regression
Learning Objectives - Session 1
By the end of today’s class, you will be able to:
- Fit simple and multiple linear regression models using
lm()in R - Interpret regression coefficients, R-squared, and p-values
- Understand the mathematical foundation of least squares estimation
- Explain properties of estimators
Our Dataset: Palmer Penguins 🐧
Today we’ll work with the famous Palmer Penguins dataset!
Research Question:
Can we describe penguin body mass using other measurements?
The Linear Regression Model
Simple Linear Regression: \[Y_i = \beta_0 + \beta_1 X_i + \epsilon_i\]
where:
- \(Y_i\) is the response variable for observation \(i\)
- \(X_i\) is the predictor variable
- \(\beta_0\) is the intercept
- \(\beta_1\) is the slope
- \(\epsilon_i\) is the error term
Multiple Linear Regression
Multiple Linear Regression: \[Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_p X_{ip} + \epsilon_i\]
Or in matrix notation: \[\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}\]
where \(\mathbf{X}\) is the design matrix and \(\boldsymbol{\beta}\) is the coefficient vector.
Assumptions of Linear Regression
Before we dive into estimation, what do we assume?
- Linearity: The relationship between \(X\) and \(Y\) is linear
- Independence: Observations are independent
- Homoscedasticity: Constant variance of errors (\(\text{Var}(\epsilon_i) = \sigma^2\))
- Normality: Errors are normally distributed (\(\epsilon_i \sim N(0, \sigma^2)\))
We write this as: \(\epsilon_i \overset{iid}{\sim} N(0, \sigma^2)\)
🎯 Exercise: Explore the Data
library(palmerpenguins)
data(penguins)
# Your turn:
# 1. What variables are available?
# 2. Which could be predictors?
# 3. Create a scatter plot of body_mass_g vs flipper_length_mm
# 4. Which other variables would you try as predictors?
ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE)Discussion: What relationships do you see?
Least Squares Estimation: The Goal
We want to find \(\hat{\beta}_0\) and \(\hat{\beta}_1\) that minimize the Sum of Squared Residuals (SSR):
\[\text{SSR} = \sum_{i=1}^{n} (Y_i - \hat{Y}_i)^2 = \sum_{i=1}^{n} (Y_i - \hat{\beta}_0 - \hat{\beta}_1 X_i)^2\]
This is called the Ordinary Least Squares (OLS) criterion.
Deriving the OLS Estimators
Take partial derivatives and set to zero:
\[\frac{\partial \text{SSR}}{\partial \beta_0} = -2\sum_{i=1}^{n}(Y_i - \hat{\beta}_0 - \hat{\beta}_1 X_i) = 0\]
\[\frac{\partial \text{SSR}}{\partial \beta_1} = -2\sum_{i=1}^{n}X_i(Y_i - \hat{\beta}_0 - \hat{\beta}_1 X_i) = 0\]
These are the normal equations.
OLS Solutions
Solving the normal equations gives:
\[\hat{\beta}_1 = \frac{\sum_{i=1}^{n}(X_i - \bar{X})(Y_i - \bar{Y})}{\sum_{i=1}^{n}(X_i - \bar{X})^2} = \frac{S_{XY}}{S_{XX}}\]
\[\hat{\beta}_0 = \bar{Y} - \hat{\beta}_1\bar{X}\]
Key insight: The regression line always passes through \((\bar{X}, \bar{Y})\)!
Matrix Form: OLS Estimator
In multiple regression, the OLS estimator is:
\[\hat{\boldsymbol{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{Y}\]
This elegant formula gives us all coefficient estimates at once!
Condition: \(\mathbf{X}^T\mathbf{X}\) must be invertible (no perfect multicollinearity).
Fitting Simple Linear Regression in R
The lm() function does all the matrix algebra for us!
🎯 Exercise: Linear Algebra
# Fit a simple linear regression model using LA
# Extract the response variable (y)
y <- penguins$body_mass_g
# Create the design matrix (X)
# Include intercept (column of 1s) and predictor variable
X <- cbind(1, penguins$flipper_length_mm)
# Remove rows with missing values (to match what lm() does)
complete_cases <- complete.cases(y, X)
y <- y[complete_cases]
X <- X[complete_cases, ]
# Calculate coefficients using linear algebra: β = (X'X)^(-1)X'y
beta_hat <- ## complete the code
# Display the results
beta_hat
# Compare with lm() results
model1 <- lm(body_mass_g ~ flipper_length_mm, data = penguins)
coef(model1)
# Give names to our manually calculated coefficients
rownames(beta_hat) <- c("(Intercept)", "flipper_length_mm")
beta_hatInterpreting the Intercept (\(\hat{\beta}_0\))
The intercept is the expected value of \(Y\) when \(X = 0\).
Example: If predicting body mass from flipper length: \[\hat{\text{Body Mass}} = -5780.83 + 49.69 \times \text{Flipper Length}\]
- \(\hat{\beta}_0 = -5780.83\) g: Expected body mass when flipper length is 0 mm
Important: This has no practical interpretation here (a penguin can’t have 0 mm flippers!).
Interpreting the Slope (\(\hat{\beta}_1\))
The slope is the expected change in \(Y\) for a one-unit increase in \(X\).
Example: \[\hat{\text{Body Mass}} = -5780.83 + 49.69 \times \text{Flipper Length}\]
- \(\hat{\beta}_1 = 49.69\) g/mm: For every 1 mm increase in flipper length, body mass increases by about 49.69 grams on average
Units matter! Always interpret in context with units.
Multiple Regression in R
Each coefficient is interpreted as the expected change in \(Y\) for a one-unit change in that predictor, holding all other predictors constant.
🎯 Exercise: Multiple Regression
Properties of OLS Estimators
Under the standard assumptions, OLS estimators are:
- Unbiased: \(E[\hat{\beta}_j] = \beta_j\)
- Consistent: \(\hat{\beta}_j \xrightarrow{p} \beta_j\) as \(n \to \infty\)
- BLUE: Best Linear Unbiased Estimator (Gauss-Markov Theorem)
BLUE means: Among all linear unbiased estimators, OLS has the smallest variance!
Variance of OLS Estimators
For simple linear regression:
\[\text{Var}(\hat{\beta}_1) = \frac{\sigma^2}{\sum_{i=1}^{n}(X_i - \bar{X})^2} = \frac{\sigma^2}{S_{XX}}\]
\[\text{Var}(\hat{\beta}_0) = \sigma^2\left(\frac{1}{n} + \frac{\bar{X}^2}{S_{XX}}\right)\]
Key insight: More spread in \(X\) (larger \(S_{XX}\)) gives more precise estimates!
Estimating \(\sigma^2\)
We estimate the error variance using the Mean Squared Error (MSE):
\[\hat{\sigma}^2 = \text{MSE} = \frac{\sum_{i=1}^{n}(Y_i - \hat{Y}_i)^2}{n - p - 1} = \frac{\text{SSR}}{n - p - 1}\]
where \(p\) is the number of predictors.
The denominator \((n - p - 1)\) gives us degrees of freedom.
R-squared: Measuring Model Fit
R-squared (\(R^2\)) measures the proportion of variance in \(Y\) explained by the model:
\[R^2 = 1 - \frac{\text{SSR}}{\text{SST}} = 1 - \frac{\sum(Y_i - \hat{Y}_i)^2}{\sum(Y_i - \bar{Y})^2}\]
- \(R^2 = 0\): Model explains nothing
- \(R^2 = 1\): Model explains everything
- Typically: \(0 < R^2 < 1\)
Adjusted R-squared
Problem: \(R^2\) always increases when we add predictors!
Solution: Adjusted R-squared penalizes model complexity:
\[R^2_{\text{adj}} = 1 - \frac{\text{SSR}/(n-p-1)}{\text{SST}/(n-1)} = 1 - \frac{n-1}{n-p-1}(1-R^2)\]
Use \(R^2_{\text{adj}}\) when comparing models with different numbers of predictors.
🎯 Exercise: Interpret R-squared
Hypothesis Testing for Coefficients
We test: \(H_0: \beta_j = 0\) vs. \(H_a: \beta_j \neq 0\)
Test statistic: \[t = \frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)} \sim t_{n-p-1}\]
p-value: Probability of observing a test statistic of this value or more extreme given that \(H_0\) is true.
Let’s talk about p-values
- In 2016, the ASA released an unprecedented statement addressing widespread misunderstanding of p-values
- P-values have contributed to a reproducibility crisis in science
- While useful, they are frequently misinterpreted and misused
Reference: Wasserstein, R. L., & Lazar, N. A. (2016). The ASA Statement on p-Values: Context, Process, and Purpose. The American Statistician, 70(2), 129-133.
What P-Values Actually Mean
ASA Principle #1
“P-values can indicate how incompatible the data are with a specified statistical model.”
A p-value tells us: The probability of obtaining results at least as extreme as observed, assuming the null hypothesis is true
A p-value does NOT tell us:
- The probability that the hypothesis is true
- The probability that results were produced by random chance alone
- The importance or magnitude of an effect
The Arbitrary Threshold Problem
ASA Principles #2 & #3
P-values do not measure the probability that the studied hypothesis is true. Scientific conclusions should not be based solely on whether a p-value passes a specific threshold.
Problems with p < 0.05:
- Dichotomania: Treating p = 0.049 as fundamentally different from p = 0.051
- Publication bias: Favoring “significant” results over null findings
- P-hacking: Manipulating analyses until significance is achieved
- Selective reporting: Hiding non-significant results
Moving Forward: Better Practice
ASA Recommendations
Proper inference requires full reporting, transparency, and consideration of multiple factors beyond p-values.
Key practices:
- Report effect sizes and confidence intervals
- Consider study design and data quality
- Account for multiple comparisons
- Replicate findings across studies
- Use domain expertise and context
- Embrace uncertainty in conclusions
Bottom line: Accept uncertainty, be thoughtful, be open, and be modest in your claims.
Confidence Intervals for Coefficients
A \((1-\alpha)100\%\) confidence interval for \(\beta_j\):
\[\hat{\beta}_j \pm t_{\alpha/2, n-p-1} \times \text{SE}(\hat{\beta}_j)\]
Interpretation: We are \((1-\alpha)100\%\) confident that the true \(\beta_j\) lies in this interval.
Getting Confidence Intervals in R
# 95% confidence intervals for coefficients
confint(model2, level = 0.95)
# For predictions
new_penguin <- data.frame(flipper_length_mm = 200,
bill_length_mm = 45,
bill_depth_mm = 15)
predict(model2, newdata = new_penguin,
interval = "confidence")
predict(model2, newdata = new_penguin,
interval = "prediction")🎯 Exercise: Hypothesis Testing
Visualization: Fitted Values
library(ggplot2)
# Add predictions to data
penguins_clean <- penguins %>%
drop_na(any_of(c("body_mass_g","bill_length_mm",
"bill_depth_mm")))%>%
mutate(fitted = fitted(model2),
residuals = residuals(model2))
# Plot actual vs fitted
ggplot(penguins_clean, aes(x = fitted, y = body_mass_g)) +
geom_point(alpha = 0.5) +
geom_abline(slope = 1, intercept = 0, color = "blue") +
labs(title = "Actual vs Fitted Values", x = "Fitted Body Mass (g)",y = "Actual Body Mass (g)")F-test for Overall Model Significance
Tests: \(H_0\): All slope coefficients are zero (\(\beta_1 = \beta_2 = \cdots = \beta_p = 0\))
Test statistic: \[F = \frac{(\text{SST} - \text{SSR})/p}{\text{SSR}/(n-p-1)} = \frac{\text{MSR}}{\text{MSE}} \sim F_{p, n-p-1}\]
This appears at the bottom of summary() output.
🎯 Exercise: Overall Model Test
Practical Tips: Standardized Coefficients
When predictors are on different scales, standardized coefficients help compare importance:
Now coefficients represent change in \(Y\) for a 1 SD change in \(X\).
Categorical Predictors
R automatically creates dummy variables for factors:
The reference category (Adelie) is absorbed into the intercept.
Interaction Terms
Test whether the effect of one predictor depends on another:
🎯 Exercise: Build Your Best Model
Session 1 Wrap-up
Key concepts covered:
- OLS estimation and its mathematical foundation
- Properties of estimators (unbiased, BLUE)
- Interpreting coefficients, R-squared, and p-values
- Simple and multiple regression in R using
lm() - Hypothesis testing for individual coefficients and overall model
Next session: Model diagnostics and selection!
Session 2: Model Diagnostics & Selection
Learning Objectives - Session 2
By the end of today’s class, you will be able to:
- Assess model assumptions using diagnostic plots
- Calculate and interpret model fit statistics
- Understand the mathematical foundation of diagnostic tests
- Perform model selection using stepwise regression, AIC, and cross-validation
- Explain the mathematical basis for model selection criteria
Why Check Assumptions?
Remember our assumptions?
- Linearity
- Independence
- Homoscedasticity (constant variance)
- Normality of errors
If assumptions are violated: Estimates may be biased, hypothesis tests invalid, predictions unreliable.
We use diagnostic plots and statistical tests to check these!
Residuals: The Key to Diagnostics
Residual for observation \(i\): \[e_i = Y_i - \hat{Y}_i\]
Standardized residual: \[r_i = \frac{e_i}{\hat{\sigma}\sqrt{1 - h_{ii}}}\]
where \(h_{ii}\) is the leverage (diagonal element of the hat matrix \(\mathbf{H} = \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\)).
Getting Residuals in R
Diagnostic Plot 1: Residuals vs Fitted
Purpose: Check linearity and homoscedasticity
Interpreting Residuals vs Fitted
What to look for:
- Random scatter around zero: Assumptions satisfied ✓
- Pattern or curve: Linearity violated ✗
- Funnel shape: Heteroscedasticity (non-constant variance) ✗
- Outliers: Points far from zero
🎯 Exercise: Residuals vs Fitted
Diagnostic Plot 2: Q-Q Plot
Purpose: Check normality of residuals
What to look for: Points should fall approximately on the diagonal line.
Interpreting Q-Q Plots
Patterns indicate:
- S-shape: Heavy tails (more extreme values than normal)
- Inverse S-shape: Light tails
- Points above line on right: Right skew
- Points below line on left: Left skew
Minor deviations: Usually okay, especially with large \(n\)
Diagnostic Plot 3: Scale-Location
Purpose: Check homoscedasticity
This plots \(\sqrt{|\text{standardized residuals}|}\) vs fitted values.
What to look for: Horizontal line with equal spread → constant variance ✓
Diagnostic Plot 4: Residuals vs Leverage
Purpose: Identify influential observations
Cook’s distance measures influence of each observation.
Rule of thumb: Cook’s distance > 1 indicates influential point.
🎯 Exercise: All Diagnostic Plots
Statistical Tests: Breusch-Pagan Test
Tests for heteroscedasticity.
\(H_0\): Homoscedasticity (constant variance)
If p-value is very small: Reject \(H_0\), evidence of heteroscedasticity.
Statistical Tests: Shapiro-Wilk Test
Tests for normality of residuals.
\(H_0\): Residuals are normally distributed
Important: With large \(n\), this test is very sensitive. Visual inspection (Q-Q plot) often more useful.
Durbin-Watson Test for Independence
Tests for autocorrelation in residuals (common in time series).
\(H_0\): No autocorrelation
Test statistic ≈ 2 suggests no autocorrelation.
🎯 Exercise: Formal Tests
What If Assumptions Are Violated?
Potential solutions:
- Transform variables (log, square root, etc.)
- Add polynomial terms (quadratic, cubic)
- Include interactions
- Remove outliers (carefully!)
- Use robust regression methods
- Try a different model family (GLM, GAM)
Model Selection: The Problem
We often have many potential predictors. Questions:
- Which predictors should we include?
- Is a more complex model better? - Description vs Prediction
- How do we balance fit vs. complexity?
Goal: Find the model that predicts best on new data!
Bias-Variance Tradeoff
Underfitting (too simple):
- High bias, low variance
- Poor fit to training and test data
Overfitting (too complex):
- Low bias, high variance
- Great fit to training data, poor on test data
Goal: Find the sweet spot!
Model Selection Strategies
- All subsets regression: Try all possible models (computationally expensive)
- Stepwise regression: Add/remove predictors sequentially
- Information criteria: AIC, BIC
- Cross-validation: Estimate prediction error directly
- Regularization: Penalize complexity (LASSO, Ridge) - next course!
Forward Selection
Algorithm:
- Start with intercept-only model
- Add predictor that improves model most
- Repeat until no improvement
- Stop based on criterion (AIC, p-value, etc.)
Backward Selection
Algorithm:
- Start with full model (all predictors)
- Remove predictor that hurts model least
- Repeat until no improvement
- Stop based on criterion
Stepwise Selection
Combines forward and backward:
- At each step, can add or remove a predictor
- More flexible than pure forward or backward
🎯 Exercise: Stepwise Selection
# Create a full model with all available predictors
full_model <- lm(body_mass_g ~ flipper_length_mm + bill_length_mm +
bill_depth_mm + species + sex + island,
data = penguins)
# Try different directions
step_forward <- step(lm(body_mass_g ~ 1, data = penguins),
scope = formula(full_model),
direction = "forward")
step_backward <- step(full_model, direction = "backward")
# Do they give the same final model?Akaike Information Criterion (AIC)
AIC balances model fit and complexity:
\[\text{AIC} = -2\log(L) + 2p\]
where \(L\) is the likelihood and \(p\) is the number of parameters.
For linear regression: \[\text{AIC} = n\log(\text{SSR}/n) + 2p\]
Lower AIC is better! AIC penalizes adding parameters.
Bayesian Information Criterion (BIC)
BIC is similar to AIC but penalizes complexity more:
\[\text{BIC} = -2\log(L) + p\log(n)\]
For linear regression: \[\text{BIC} = n\log(\text{SSR}/n) + p\log(n)\]
Lower BIC is better! BIC tends to select simpler models than AIC.
Comparing Models with AIC/BIC
# Fit several models
model1 <- lm(body_mass_g ~ flipper_length_mm, data = penguins)
model2 <- lm(body_mass_g ~ flipper_length_mm + bill_length_mm, data = penguins)
model3 <- lm(body_mass_g ~ flipper_length_mm + bill_length_mm + species,
data = penguins)
# Compare AIC
AIC(model1, model2, model3)
# Compare BIC
BIC(model1, model2, model3)Guideline: Difference in AIC > 2 suggests models are meaningfully different.
🎯 Exercise: Compare with AIC/BIC
Cross-Validation: The Gold Standard
Idea: Estimate prediction error on unseen data.
K-Fold Cross-Validation:
- Split data into \(K\) folds
- For each fold: train on \(K-1\) folds, test on remaining fold
- Average the prediction errors
Common: \(K = 5\) or \(K = 10\)
Leave-One-Out Cross-Validation (LOOCV)
Special case: \(K = n\)
CV error: \[\text{CV} = \frac{1}{n}\sum_{i=1}^{n}(Y_i - \hat{Y}_{-i})^2\]
where \(\hat{Y}_{-i}\) is prediction for \(i\) using model trained without observation \(i\).
Computationally expensive but can use shortcut for linear regression!
LOOCV Shortcut for Linear Regression
\[\text{CV} = \frac{1}{n}\sum_{i=1}^{n}\left(\frac{e_i}{1 - h_{ii}}\right)^2\] where \(e_i\) is the residual and \(h_{ii}\) is the leverage.
This is called the PRESS statistic (Predicted Residual Sum of Squares).
Best Practices for Model Selection
- Start with theory: Use domain knowledge
- Check assumptions: Before comparing models
- Use multiple criteria: Don’t rely on just one metric
- Validate: Always test on new data
- Keep it simple: Prefer simpler models when performance is similar
- Document: Record your selection process
Common Pitfalls to Avoid
❌ Don’t:
- Test on the same data used for selection
- Use automated selection blindly
- Ignore assumption violations
- Focus only on R-squared
- Include perfectly correlated predictors
✓ Do:
- Split data or use cross-validation
- Understand why predictors are selected
- Check diagnostics for final model
- Consider multiple criteria
- Check for multicollinearity
🎯 Final Exercise: Complete Analysis
Conduct a complete analysis:
- Fit several candidate models
- Check assumptions with diagnostic plots
- Calculate VIF for multicollinearity
- Compare models using AIC, BIC
- Perform cross-validation
- Select your final model
- Interpret coefficients in context
- Make predictions for new penguins
Resources for Further Learning
📚 Recommended:
- James et al. (2021) - An Introduction to Statistical Learning
- Faraway (2014) - Linear Models with R
💻 R packages:
- car: Diagnostics and VIF
- lmtest: Statistical tests
- caret: Cross-validation
- MASS: Stepwise selection
STAT 204 – Intro to Statistical Data Analysis