Regression analysis is one of the most fundamental and widely used statistical techniques in data science, machine learning, and scientific research. At its core, regression helps us understand relationships between variables and make predictions based on those relationships.
Formal Definition: Regression is a statistical method for estimating the relationships between a dependent variable (often called the response or outcome) and one or more independent variables (called predictors, features, or explanatory variables).
The term “regression” was coined by Sir Francis Galton in the 19th century. While studying the heights of parents and their children, Galton observed that children’s heights tended to “regress” toward the average height of the population. This phenomenon, called regression toward the mean, gave birth to the statistical technique we use today.
Imagine you are a botanist studying iris flowers. You might ask:
Regression helps answer these questions quantitatively. It allows us to:
Before diving into regression, let’s understand our dataset. The Iris dataset is perhaps the most famous dataset in statistics and machine learning, introduced by Sir Ronald A. Fisher in his 1936 paper on discriminant analysis.
The dataset contains measurements of 150 iris flowers from three species: - Iris setosa (50 flowers) - Iris versicolor (50 flowers) - Iris virginica (50 flowers)
Understanding flower anatomy is crucial for understanding this dataset:
| Term | Description | Function |
|---|---|---|
| Sepal | The outer, leaf-like protective structures of a flower bud | Protect the flower before it blooms; often green and photosynthetic |
| Petal | The inner, often colorful structures of a flower | Attract pollinators with color and scent |
In the iris flower: - Sepals are the three outer segments that hang downward (called “falls” in iris terminology) - Petals are the three inner segments that stand upright (called “standards”)
Figure 1: Anatomy of an iris flower showing petals (standards) and sepals (falls). Source: USDA Forest Service.
Figure 2: The three iris species in our dataset. Notice the differences in petal and sepal sizes. Source: GitHub educational repository.
For each flower, four measurements were recorded (in centimeters):
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
## 7 4.6 3.4 1.4 0.3 setosa
## 8 5.0 3.4 1.5 0.2 setosa
## 9 4.4 2.9 1.4 0.2 setosa
## 10 4.9 3.1 1.5 0.1 setosa
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
##
## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
Key Observation: The iris dataset has 150 observations and 5 variables – 4 numeric measurements and 1 categorical variable (Species) with 3 levels.
# Count flowers per species
ggplot(iris, aes(x = Species, fill = Species)) +
geom_bar(color = "white", linewidth = 1.5) +
scale_fill_manual(values = iris_colors) +
geom_text(
stat = "count", aes(label = after_stat(count)),
vjust = -0.5, size = 6, fontface = "bold"
) +
labs(
title = "Distribution of Iris Species",
subtitle = "50 flowers from each of the three species",
x = "Species",
y = "Count"
) +
theme(legend.position = "none") +
ylim(0, 60)# Create a scatterplot matrix to see all relationships
pairs(iris[, 1:4],
col = iris_colors[iris$Species],
pch = 19,
main = "Scatterplot Matrix of Iris Measurements",
cex.labels = 1.5
)Simple Linear Regression (SLR) is the simplest form of regression. It models the relationship between two continuous variables using a straight line.
\[Y = \beta_0 + \beta_1 X + \varepsilon\]
Where: - \(Y\) = Dependent variable (response) - \(X\) = Independent variable (predictor) - \(\beta_0\) = Intercept (value of Y when X = 0) - \(\beta_1\) = Slope (change in Y for a one-unit change in X) - \(\varepsilon\) = Error term (random noise)
The Goal: Find the “best-fit” line that minimizes the sum of squared errors (differences between observed and predicted values). This is called the Least Squares Method.
Let’s predict Petal.Length (Y) using Sepal.Length (X).
# Fit a simple linear regression model
model_slr <- lm(Petal.Length ~ Sepal.Length, data = iris)
# View the model summary
summary(model_slr)##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.47747 -0.59072 -0.00668 0.60484 2.49512
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
## Sepal.Length 1.85843 0.08586 21.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8678 on 148 degrees of freedom
## Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
## F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16
The output shows several important pieces of information:
# Visualize the regression line
ggplot(iris, aes(x = Sepal.Length, y = Petal.Length)) +
geom_point(aes(color = Species), size = 3, alpha = 0.8) +
geom_smooth(
method = "lm", se = TRUE, color = "black",
linetype = "dashed", linewidth = 1.5
) +
scale_color_manual(values = iris_colors) +
labs(
title = "Simple Linear Regression: Petal Length vs Sepal Length",
subtitle = "Black dashed line shows the best-fit regression line",
x = "Sepal Length (cm)",
y = "Petal Length (cm)"
)Interpretation: For every 1 cm increase in Sepal Length, Petal Length increases by approximately 1.86 cm (the slope coefficient). The R-squared value tells us what proportion of variation in Petal Length is explained by Sepal Length.
Multiple Linear Regression (MLR) extends simple linear regression by using two or more predictors to model the response variable.
\[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p + \varepsilon\]
Where: - \(X_1, X_2, \ldots, X_p\) are the predictor variables - \(\beta_1, \beta_2, \ldots, \beta_p\) are the partial regression coefficients
# Fit multiple linear regression
model_mlr <- lm(Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width,
data = iris
)
# View summary
summary(model_mlr)##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width,
## data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99333 -0.17656 -0.01004 0.18558 1.06909
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.26271 0.29741 -0.883 0.379
## Sepal.Length 0.72914 0.05832 12.502 <2e-16 ***
## Sepal.Width -0.64601 0.06850 -9.431 <2e-16 ***
## Petal.Width 1.44679 0.06761 21.399 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.319 on 146 degrees of freedom
## Multiple R-squared: 0.968, Adjusted R-squared: 0.9674
## F-statistic: 1473 on 3 and 146 DF, p-value: < 2.2e-16
In multiple regression, each coefficient represents the change in Y for a one-unit change in that predictor, holding all other predictors constant. This is why they’re called partial regression coefficients.
# Visualize actual vs predicted values
iris$predicted_mlr <- predict(model_mlr)
ggplot(iris, aes(x = predicted_mlr, y = Petal.Length)) +
geom_point(aes(color = Species), size = 3, alpha = 0.8) +
geom_abline(
intercept = 0, slope = 1, color = "black",
linetype = "dashed", linewidth = 1.5
) +
scale_color_manual(values = iris_colors) +
labs(
title = "Multiple Regression: Actual vs Predicted Petal Length",
subtitle = "Points on the dashed line indicate perfect predictions",
x = "Predicted Petal Length (cm)",
y = "Actual Petal Length (cm)"
)# Compare R-squared values
cat(
"Simple Linear Regression R-squared:",
round(summary(model_slr)$r.squared, 4), "\n"
)## Simple Linear Regression R-squared: 0.76
## Multiple Linear Regression R-squared: 0.968
Key Insight: Multiple regression typically has higher R-squared because adding more predictors always improves fit. However, we should use Adjusted R-squared for fair comparison, as it penalizes unnecessary predictors.
Sometimes relationships are not linear. Polynomial regression models curved relationships by including polynomial terms.
\[Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \beta_3 X^3 + \cdots + \varepsilon\]
# Fit polynomial regression (degree 2)
model_poly <- lm(Petal.Length ~ Sepal.Length + I(Sepal.Length^2), data = iris)
# Compare with linear model
summary(model_poly)##
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + I(Sepal.Length^2),
## data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6751 -0.5138 0.1218 0.5356 2.6287
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.44671 3.08770 -5.650 8.04e-08 ***
## Sepal.Length 5.39216 1.04465 5.162 7.81e-07 ***
## I(Sepal.Length^2) -0.29586 0.08719 -3.393 0.000887 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8385 on 147 degrees of freedom
## Multiple R-squared: 0.7774, Adjusted R-squared: 0.7744
## F-statistic: 256.7 on 2 and 147 DF, p-value: < 2.2e-16
# Visualize polynomial fit
ggplot(iris, aes(x = Sepal.Length, y = Petal.Length)) +
geom_point(aes(color = Species), size = 3, alpha = 0.8) +
geom_smooth(
method = "lm", formula = y ~ poly(x, 2),
color = "black", se = TRUE, linewidth = 1.5
) +
scale_color_manual(values = iris_colors) +
labs(
title = "Polynomial Regression (Degree 2)",
subtitle = "Curved line captures non-linear relationships",
x = "Sepal Length (cm)",
y = "Petal Length (cm)"
)Logistic regression is used when the response variable is categorical (usually binary). Instead of predicting a continuous value, we predict the probability of belonging to a category.
\[P(Y = 1) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p)}}\]
The logistic function ensures predicted probabilities are between 0 and 1.
# Create binary response: Is it setosa?
iris$is_setosa <- ifelse(iris$Species == "setosa", 1, 0)
# Fit logistic regression
model_logit <- glm(
is_setosa ~ Sepal.Length + Sepal.Width +
Petal.Length + Petal.Width,
data = iris, family = binomial
)
# View summary
summary(model_logit)##
## Call:
## glm(formula = is_setosa ~ Sepal.Length + Sepal.Width + Petal.Length +
## Petal.Width, family = binomial, data = iris)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -16.946 457457.097 0 1
## Sepal.Length 11.759 130504.042 0 1
## Sepal.Width 7.842 59415.385 0 1
## Petal.Length -20.088 107724.594 0 1
## Petal.Width -21.608 154350.616 0 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1.9095e+02 on 149 degrees of freedom
## Residual deviance: 3.2940e-09 on 145 degrees of freedom
## AIC: 10
##
## Number of Fisher Scoring iterations: 25
# Predict probabilities
iris$prob_setosa <- predict(model_logit, type = "response")
# Visualize
ggplot(iris, aes(x = Petal.Length, y = prob_setosa)) +
geom_point(aes(color = Species), size = 3, alpha = 0.8) +
geom_smooth(
method = "glm", method.args = list(family = binomial),
color = "black", se = TRUE, linewidth = 1.5
) +
scale_color_manual(values = iris_colors) +
labs(
title = "Logistic Regression: Probability of Being Setosa",
subtitle = "S-shaped curve shows probability transition",
x = "Petal Length (cm)",
y = "Probability of Setosa"
)Linear Discriminant Analysis (LDA) is both a dimensionality reduction technique and a classification method. Developed by Fisher in 1936 (yes, the same Fisher who gave us the iris dataset!), LDA finds linear combinations of features that best separate classes.
LDA makes three key assumptions: 1. Normality: Each class follows a multivariate normal distribution 2. Equal covariance matrices: All classes share the same covariance structure 3. Independence: Observations are independent
The algorithm works in two steps:
LDA seeks to find a projection vector w that maximizes:
\[J(w) = \frac{w^T S_B w}{w^T S_W w}\]
Where: - \(S_B\) = Between-class scatter matrix (measures separation between classes) - \(S_W\) = Within-class scatter matrix (measures spread within each class)
Intuition: LDA finds the “best viewing angle” to separate the species. Imagine rotating a 3D scatter plot until the three species are most clearly separated – that’s what LDA does mathematically.
# Perform LDA
model_lda <- lda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, data = iris)
# View the model
model_lda## Call:
## lda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
## data = iris)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 5.006 3.428 1.462 0.246
## versicolor 5.936 2.770 4.260 1.326
## virginica 6.588 2.974 5.552 2.026
##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.8293776 -0.02410215
## Sepal.Width 1.5344731 -2.16452123
## Petal.Length -2.2012117 0.93192121
## Petal.Width -2.8104603 -2.83918785
##
## Proportion of trace:
## LD1 LD2
## 0.9912 0.0088
The output shows: - Prior probabilities: Proportion of each species in the data - Group means: Average measurements for each species - Coefficients of linear discriminants: Weights for each variable in creating the discriminant axes
# Predict LDA scores
lda_scores <- predict(model_lda)$x
# Create a data frame for plotting
lda_df <- data.frame(
LD1 = lda_scores[, 1],
LD2 = lda_scores[, 2],
Species = iris$Species
)
# Plot LDA projection
ggplot(lda_df, aes(x = LD1, y = LD2, color = Species)) +
geom_point(size = 4, alpha = 0.8) +
stat_ellipse(level = 0.95, linewidth = 1.5) +
scale_color_manual(values = iris_colors) +
labs(
title = "LDA Projection of Iris Species",
subtitle = "Two linear discriminants separate the three species",
x = "Linear Discriminant 1 (LD1)",
y = "Linear Discriminant 2 (LD2)"
) +
theme(legend.position = "right")Amazing Result: Notice how well-separated the species are in LDA space! Even though we started with 4 dimensions, just 2 discriminant axes capture most of the class separation.
# Predict species
lda_pred <- predict(model_lda)$class
# Confusion matrix
confusion_lda <- table(Predicted = lda_pred, Actual = iris$Species)
print(confusion_lda)## Actual
## Predicted setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 48 1
## virginica 0 2 49
# Calculate accuracy
accuracy_lda <- sum(diag(confusion_lda)) / sum(confusion_lda)
cat("\nLDA Classification Accuracy:", round(accuracy_lda * 100, 2), "%\n")##
## LDA Classification Accuracy: 98 %
# Convert confusion matrix to data frame
confusion_df <- as.data.frame(confusion_lda)
ggplot(confusion_df, aes(x = Actual, y = Predicted, fill = Freq)) +
geom_tile(color = "white", linewidth = 2) +
geom_text(aes(label = Freq), size = 8, fontface = "bold", color = "white") +
scale_fill_gradient(low = "#3498DB", high = "#E74C3C") +
labs(
title = "LDA Confusion Matrix",
subtitle = paste("Accuracy:", round(accuracy_lda * 100, 2), "%"),
x = "Actual Species",
y = "Predicted Species"
) +
theme(legend.position = "none")Ridge regression adds a penalty term to prevent overfitting when predictors are highly correlated.
\[\text{Minimize: } \sum_{i=1}^{n}(y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} \beta_j^2\]
# Ridge regression using glmnet
library(glmnet)
# Prepare data
x <- as.matrix(iris[, 1:4])
y <- iris$Petal.Length
# Fit ridge regression
model_ridge <- glmnet(x, y, alpha = 0)
# Plot coefficient paths
plot(model_ridge, xvar = "lambda", label = TRUE)
title("Ridge Regression: Coefficient Paths", line = 2.5)Lasso (Least Absolute Shrinkage and Selection Operator) can shrink some coefficients to exactly zero, performing variable selection.
\[\text{Minimize: } \sum_{i=1}^{n}(y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} |\beta_j|\]
# Fit lasso regression
model_lasso <- glmnet(x, y, alpha = 1)
# Plot coefficient paths
plot(model_lasso, xvar = "lambda", label = TRUE)
title("Lasso Regression: Coefficient Paths", line = 2.5)Quantile regression estimates the conditional median or other quantiles of the response variable, rather than the mean.
library(quantreg)
# Fit quantile regression for median (0.5)
model_quant <- rq(Petal.Length ~ Sepal.Length, data = iris, tau = 0.5)
# Compare with linear regression
summary(model_quant)##
## Call: rq(formula = Petal.Length ~ Sepal.Length, tau = 0.5, data = iris)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) -6.86364 -9.36601 -6.24413
## Sepal.Length 1.81818 1.70554 2.13286
library(ggplot2)
library(quantreg)
# Visualize multiple quantile regressions with distinct colors
ggplot(iris, aes(x = Sepal.Length, y = Petal.Length)) +
geom_point(aes(color = Species), size = 3, alpha = 0.7) +
scale_color_manual(values = iris_colors) +
# Use ggnewscale or map quantiles inside aes() to separate scale mappings
geom_quantile(
aes(color = factor(after_stat(quantile))),
quantiles = c(0.1, 0.5, 0.9),
linewidth = 1.5
) +
scale_color_manual(
values = c("0.1" = "blue", "0.5" = "black", "0.9" = "red"),
labels = c("10th Percentile", "50th Percentile (Median)", "90th Percentile"),
name = "Quantile"
) +
labs(
title = "Quantile Regression at Different Levels",
subtitle = "Blue: 10th percentile, Black: median, Red: 90th percentile",
x = "Sepal Length (cm)",
y = "Petal Length (cm)"
)| Method | Type | Response Variable | Key Feature | Best Use Case |
|---|---|---|---|---|
| Simple Linear | Regression | Continuous | One predictor | Basic relationships |
| Multiple Linear | Regression | Continuous | Multiple predictors | Complex relationships |
| Polynomial | Regression | Continuous | Curved relationships | Non-linear patterns |
| Logistic | Classification | Binary/Categorical | Probability estimation | Binary outcomes |
| LDA | Classification | Categorical | Dimension reduction | Multi-class problems |
| Ridge | Regression | Continuous | L2 regularization | Correlated predictors |
| Lasso | Regression | Continuous | L1 regularization | Feature selection |
| Quantile | Regression | Continuous | Robust to outliers | Non-normal errors |
# Create a comparison of R-squared values
methods <- c("Simple Linear", "Multiple Linear", "Polynomial")
r_squared <- c(
summary(model_slr)$r.squared,
summary(model_mlr)$r.squared,
summary(model_poly)$r.squared
)
comparison_df <- data.frame(Method = methods, R_squared = r_squared)
ggplot(comparison_df, aes(x = reorder(Method, R_squared), y = R_squared, fill = R_squared)) +
geom_col(color = "white", linewidth = 1.5) +
geom_text(aes(label = round(R_squared, 4)), hjust = -0.2, size = 6, fontface = "bold") +
scale_fill_gradient(low = "#3498DB", high = "#E74C3C") +
coord_flip() +
labs(
title = "Model Comparison: R-squared Values",
subtitle = "Higher values indicate better fit",
x = "Method",
y = "R-squared"
) +
theme(legend.position = "none") +
ylim(0, 1.1)Common Mistake: Students often use complex methods when simple ones would suffice. Always start with the simplest method and add complexity only when necessary.
# Function to check assumptions
check_assumptions <- function(model) {
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))
}
# Example usage
check_assumptions(model_mlr)Remember: Statistical significance does not imply practical significance. Always consider effect sizes and confidence intervals.
Regression is about relationships: It helps us understand how variables relate to each other and make predictions.
Start simple: Begin with simple linear regression and add complexity only when needed.
Check assumptions: Always verify that your data meets the assumptions of your chosen method.
LDA is powerful: For classification problems with continuous predictors, LDA often works remarkably well (as demonstrated by our 98% accuracy on iris data!).
No single best method: The “best” method depends on your data, your question, and your constraints.
The iris dataset has taught generations of statisticians and data scientists because it: - Is small enough to understand completely - Has clear structure and patterns - Demonstrates both regression and classification - Shows the power of dimensionality reduction (LDA) - Connects to real-world biology
# ============================================================
# COMPLETE CODE FOR ALL METHODS
# ============================================================
# 1. Load libraries
library(ggplot2)
library(dplyr)
library(MASS)
library(glmnet)
library(quantreg)
# 2. Load data
data(iris)
# 3. Simple Linear Regression
model_slr <- lm(Petal.Length ~ Sepal.Length, data = iris)
summary(model_slr)
# 4. Multiple Linear Regression
model_mlr <- lm(Petal.Length ~ Sepal.Length + Sepal.Width + Petal.Width,
data = iris
)
summary(model_mlr)
# 5. Polynomial Regression
model_poly <- lm(Petal.Length ~ Sepal.Length + I(Sepal.Length^2), data = iris)
summary(model_poly)
# 6. Logistic Regression
iris$is_setosa <- ifelse(iris$Species == "setosa", 1, 0)
model_logit <- glm(
is_setosa ~ Sepal.Length + Sepal.Width +
Petal.Length + Petal.Width,
data = iris, family = binomial
)
summary(model_logit)
# 7. LDA
model_lda <- lda(Species ~ ., data = iris)
predictions <- predict(model_lda)$class
accuracy <- mean(predictions == iris$Species)
print(paste("LDA Accuracy:", round(accuracy * 100, 2), "%"))
# 8. Ridge Regression
x <- as.matrix(iris[, 1:4])
y <- iris$Petal.Length
model_ridge <- glmnet(x, y, alpha = 0)
# 9. Lasso Regression
model_lasso <- glmnet(x, y, alpha = 1)
# 10. Quantile Regression
model_quant <- rq(Petal.Length ~ Sepal.Length, data = iris, tau = 0.5)
summary(model_quant)Document prepared for undergraduate teaching in Civil Engineering and Data Science courses.