Linear Regression
Model linear relationships between input features and a continuous target using scikit-learn.
Concept Overview
Linear regression is one of the most fundamental algorithms in machine learning. It fits a line (or hyperplane for multiple features) that minimizes the sum of squared residuals between predictions and actual values.
For simple linear regression with one feature, the model is:
ŷ = β₀ + β₁·x
For multiple linear regression with n features, the model is:
ŷ = β₀ + β₁x₁ + β₂x₂ + … + βₙxₙ
Where:
- ŷ (y-hat) is the predicted value
- β₀ is the y-intercept (bias term)
- β₁, β₂, …, βₙ are the coefficients (slopes) for each feature
- x₁, x₂, …, xₙ are the input features
The goal is to find the coefficients that minimize the cost function (Mean Squared Error):
MSE = (1/n) · Σᵢ (yᵢ − ŷᵢ)²
Assumptions: Linear regression assumes:
- Linearity: The relationship between features and target is linear
- Independence: Observations are independent
- Homoscedasticity: Errors have constant variance
- Normality: Errors are normally distributed (for inference)
- No multicollinearity: Features are not highly correlated
Simple Linear Regression
Simple linear regression models the relationship between a single feature and the target variable. Let's start with a basic example.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.model_selection import train_test_split
# Generate sample data
np.random.seed(42)
sqft = np.random.normal(1500, 500, 100)
price = 100 * sqft + np.random.normal(0, 50000, 100)
df = pd.DataFrame({'sqft': sqft, 'price': price})
# Prepare data
X = df[['sqft']]
y = df['price']
# Split into train and test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Create and train the model
model = LinearRegression()
model.fit(X_train, y_train)
# Make predictions
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)
# Model parameters
print(f"Intercept (β₀): {model.intercept_:.2f}")
print(f"Coefficient (β₁): {model.coef_[0]:.2f}")
print(f"\nModel equation: price = {model.intercept_:.2f} + {model.coef_[0]:.2f} × sqft")
# Evaluate metrics
train_mse = mean_squared_error(y_train, y_train_pred)
test_mse = mean_squared_error(y_test, y_test_pred)
train_r2 = r2_score(y_train, y_train_pred)
test_r2 = r2_score(y_test, y_test_pred)
train_mae = mean_absolute_error(y_train, y_train_pred)
test_mae = mean_absolute_error(y_test, y_test_pred)
print(f"\nTraining Metrics:")
print(f" MSE: {train_mse:.2f}")
print(f" MAE: {train_mae:.2f}")
print(f" R²: {train_r2:.4f}")
print(f"\nTest Metrics:")
print(f" MSE: {test_mse:.2f}")
print(f" MAE: {test_mae:.2f}")
print(f" R²: {test_r2:.4f}")
# Visualize
plt.figure(figsize=(12, 5))
# Plot 1: Training data
plt.subplot(1, 2, 1)
plt.scatter(X_train, y_train, alpha=0.6, s=50, label='Training data')
plt.plot(X_train, y_train_pred, color='red', linewidth=2, label='Regression line')
plt.xlabel('Square Feet', fontsize=12)
plt.ylabel('Price ($)', fontsize=12)
plt.title(f'Training Set (R² = {train_r2:.4f})', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
# Plot 2: Test data
plt.subplot(1, 2, 2)
plt.scatter(X_test, y_test, alpha=0.6, s=50, label='Test data', color='green')
plt.plot(X_test, y_test_pred, color='red', linewidth=2, label='Regression line')
plt.xlabel('Square Feet', fontsize=12)
plt.ylabel('Price ($)', fontsize=12)
plt.title(f'Test Set (R² = {test_r2:.4f})', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()Always split data into train/test sets before evaluating to avoid overly optimistic metrics. The test set should only be used for final evaluation, not for model selection or hyperparameter tuning.
Multiple Linear Regression
Multiple linear regression extends simple linear regression to include multiple features, allowing us to model more complex relationships.
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split
# Generate sample data with multiple features
np.random.seed(42)
n = 200
sqft = np.random.normal(1500, 500, n)
bedrooms = np.random.randint(1, 6, n)
age = np.random.normal(10, 5, n)
price = 100 * sqft + 20000 * bedrooms - 5000 * age + np.random.normal(0, 50000, n)
df = pd.DataFrame({
'sqft': sqft,
'bedrooms': bedrooms,
'age': age,
'price': price
})
# Prepare data
X = df[['sqft', 'bedrooms', 'age']]
y = df['price']
# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Create and train model
model = LinearRegression()
model.fit(X_train, y_train)
# Make predictions
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)
# Model parameters
print("Model Coefficients:")
print(f"Intercept (β₀): {model.intercept_:.2f}")
for i, feature in enumerate(X.columns):
print(f"{feature} (β{i+1}): {model.coef_[i]:.2f}")
print(f"\nModel equation:")
print(f"price = {model.intercept_:.2f}", end="")
for i, feature in enumerate(X.columns):
sign = "+" if model.coef_[i] >= 0 else ""
print(f" {sign}{model.coef_[i]:.2f} × {feature}", end="")
print()
# Evaluate
train_r2 = r2_score(y_train, y_train_pred)
test_r2 = r2_score(y_test, y_test_pred)
train_mse = mean_squared_error(y_train, y_train_pred)
test_mse = mean_squared_error(y_test, y_test_pred)
print(f"\nTraining R²: {train_r2:.4f}")
print(f"Test R²: {test_r2:.4f}")
print(f"Training MSE: {train_mse:.2f}")
print(f"Test MSE: {test_mse:.2f}")
# Feature importance (using coefficient magnitudes when features are standardized)
print("\nFeature Importance (coefficient magnitudes):")
feature_importance = pd.DataFrame({
'Feature': X.columns,
'Coefficient': model.coef_,
'Abs_Coefficient': np.abs(model.coef_)
}).sort_values('Abs_Coefficient', ascending=False)
print(feature_importance)Multiple linear regression allows modeling relationships with multiple features simultaneously. Note that coefficient magnitudes can be misleading if features have different scales—always consider feature scaling when comparing importance.
Interpreting Coefficients
Coefficients in linear regression represent the expected change in the target variable (y) for a one-unit change in the feature, while holding all other features constant. Understanding coefficients is crucial for model interpretation.
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
# Generate sample data
np.random.seed(42)
sqft = np.random.normal(1500, 500, 100)
bedrooms = np.random.randint(1, 6, 100)
price = 100 * sqft + 20000 * bedrooms + np.random.normal(0, 50000, 100)
df = pd.DataFrame({
'sqft': sqft,
'bedrooms': bedrooms,
'price': price
})
# Model with original scale
X_original = df[['sqft', 'bedrooms']]
y = df['price']
model_original = LinearRegression()
model_original.fit(X_original, y)
print("Model with Original Scale:")
print(f"Intercept: {model_original.intercept_:.2f}")
print(f"sqft coefficient: {model_original.coef_[0]:.2f}")
print(f"bedrooms coefficient: {model_original.coef_[1]:.2f}")
print(f"\nInterpretation:")
print(f"- For each additional square foot, price increases by ${model_original.coef_[0]:.2f}")
print(f"- For each additional bedroom, price increases by ${model_original.coef_[1]:.2f}")
# Model with standardized features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_original)
model_scaled = LinearRegression()
model_scaled.fit(X_scaled, y)
print(f"\n\nModel with Standardized Features:")
print(f"Intercept: {model_scaled.intercept_:.2f}")
print(f"sqft coefficient (standardized): {model_scaled.coef_[0]:.2f}")
print(f"bedrooms coefficient (standardized): {model_scaled.coef_[1]:.2f}")
print(f"\nInterpretation (standardized):")
print(f"- A 1 standard deviation increase in sqft changes price by {model_scaled.coef_[0]:.2f} standard deviations")
print(f"- A 1 standard deviation increase in bedrooms changes price by {model_scaled.coef_[1]:.2f} standard deviations")
# Compare feature importance
print(f"\n\nFeature Importance (coefficient magnitude):")
print(f"Original scale - sqft: {abs(model_original.coef_[0]):.2f}, bedrooms: {abs(model_original.coef_[1]):.2f}")
print(f"Standardized scale - sqft: {abs(model_scaled.coef_[0]):.2f}, bedrooms: {abs(model_scaled.coef_[1]):.2f}")
# Coefficient statistics (using statsmodels for more detailed info)
try:
import statsmodels.api as sm
X_with_const = sm.add_constant(X_original)
model_stats = sm.OLS(y, X_with_const).fit()
print("\n\nDetailed Coefficient Statistics:")
print(model_stats.summary().tables[1])
except ImportError:
print("\nNote: Install statsmodels for detailed coefficient statistics (p-values, confidence intervals)")Key Points:
- Original Scale: Coefficients show the change in target per unit change in feature.
- Standardized Scale: Coefficients are comparable and show importance relative to feature variance.
- Feature Scaling: When features have different scales, coefficients can't be directly compared for importance.
- Ceteris Paribus: Coefficients assume all other features are held constant.
- Significance: Check p-values to determine if coefficients are statistically significant.
Always examine feature scaling and correlation before interpreting coefficients to avoid misleading conclusions.
Residual Analysis
Residual analysis is crucial for diagnosing model assumptions. Residuals are the differences between actual and predicted values: eᵢ = yᵢ − ŷᵢ.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from scipy import stats
# Generate sample data
np.random.seed(42)
sqft = np.random.normal(1500, 500, 100)
price = 100 * sqft + np.random.normal(0, 50000, 100)
df = pd.DataFrame({'sqft': sqft, 'price': price})
X = df[['sqft']]
y = df['price']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Train model
model = LinearRegression()
model.fit(X_train, y_train)
# Make predictions
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)
# Calculate residuals
train_residuals = y_train - y_train_pred
test_residuals = y_test - y_test_pred
# Residual plots
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
# 1. Residuals vs Predicted Values
axes[0, 0].scatter(y_train_pred, train_residuals, alpha=0.6, s=50)
axes[0, 0].axhline(y=0, color='red', linestyle='--', linewidth=2)
axes[0, 0].set_xlabel('Predicted Values', fontsize=12)
axes[0, 0].set_ylabel('Residuals', fontsize=12)
axes[0, 0].set_title('Residuals vs Predicted (Training)', fontsize=14)
axes[0, 0].grid(True, alpha=0.3)
# Random scatter indicates good fit, patterns suggest non-linearity or heteroscedasticity
# 2. Q-Q Plot (check normality)
from scipy import stats
stats.probplot(train_residuals, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot (Check Normality)', fontsize=14)
axes[0, 1].grid(True, alpha=0.3)
# Points should follow the diagonal line if residuals are normally distributed
# 3. Residuals Histogram
axes[1, 0].hist(train_residuals, bins=20, edgecolor='black', alpha=0.7)
axes[1, 0].set_xlabel('Residuals', fontsize=12)
axes[1, 0].set_ylabel('Frequency', fontsize=12)
axes[1, 0].set_title('Residuals Distribution', fontsize=14)
axes[1, 0].axvline(x=0, color='red', linestyle='--', linewidth=2)
axes[1, 0].grid(True, alpha=0.3)
# Should be approximately normally distributed with mean near 0
# 4. Residuals vs Features
axes[1, 1].scatter(X_train['sqft'], train_residuals, alpha=0.6, s=50)
axes[1, 1].axhline(y=0, color='red', linestyle='--', linewidth=2)
axes[1, 1].set_xlabel('Square Feet', fontsize=12)
axes[1, 1].set_ylabel('Residuals', fontsize=12)
axes[1, 1].set_title('Residuals vs Feature', fontsize=14)
axes[1, 1].grid(True, alpha=0.3)
# Random scatter indicates linear relationship is appropriate
plt.tight_layout()
plt.show()
# Statistical tests
print("Residual Analysis:")
print(f"Mean of residuals: {train_residuals.mean():.4f} (should be ~0)")
print(f"Standard deviation: {train_residuals.std():.4f}")
# Test for normality (Shapiro-Wilk test)
from scipy.stats import shapiro, normaltest
if len(train_residuals) <= 5000:
stat, p_value = shapiro(train_residuals)
print(f"\nShapiro-Wilk normality test:")
print(f" Statistic: {stat:.4f}")
print(f" P-value: {p_value:.6f}")
print(f" Residuals are normal: {p_value > 0.05}")
# D'Agostino-Pearson test
stat_dp, p_value_dp = normaltest(train_residuals)
print(f"\nD'Agostino-Pearson normality test:")
print(f" Statistic: {stat_dp:.4f}")
print(f" P-value: {p_value_dp:.6f}")
print(f" Residuals are normal: {p_value_dp > 0.05}")
# Check for heteroscedasticity (non-constant variance)
# Breusch-Pagan test (requires statsmodels)
try:
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan
X_with_const = sm.add_constant(X_train)
model_sm = sm.OLS(y_train, X_with_const).fit()
_, p_value_bp, _, _ = het_breuschpagan(model_sm.resid, X_with_const)
print(f"\nBreusch-Pagan test for heteroscedasticity:")
print(f" P-value: {p_value_bp:.6f}")
print(f" Constant variance (homoscedasticity): {p_value_bp > 0.05}")
except ImportError:
print("\nNote: Install statsmodels for heteroscedasticity tests")What to Look For:
- Residuals vs Predicted: Random scatter = good fit; patterns indicate non-linearity or heteroscedasticity
- Q-Q Plot: Points following diagonal = normality; deviations suggest non-normal errors
- Histogram: Bell-shaped distribution centered at 0 = normality
- Residuals vs Features: Random scatter = linear relationship is appropriate
Multicollinearity Detection
Multicollinearity occurs when features are highly correlated, making coefficient estimates unstable and difficult to interpret. Variance Inflation Factor (VIF) is a common measure.
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.api as sm
# Generate data with multicollinearity
np.random.seed(42)
sqft = np.random.normal(1500, 500, 100)
# Create correlated features
bedrooms = sqft / 300 + np.random.normal(0, 0.5, 100) # Highly correlated with sqft
bathrooms = bedrooms * 0.7 + np.random.normal(0, 0.3, 100) # Correlated with bedrooms
age = np.random.normal(10, 5, 100) # Independent
price = 100 * sqft + 20000 * bedrooms - 5000 * age + np.random.normal(0, 50000, 100)
df = pd.DataFrame({
'sqft': sqft,
'bedrooms': bedrooms,
'bathrooms': bathrooms,
'age': age,
'price': price
})
# Correlation matrix
print("Correlation Matrix:")
correlation_matrix = df[['sqft', 'bedrooms', 'bathrooms', 'age']].corr()
print(correlation_matrix)
# Visualize correlation
import matplotlib.pyplot as plt
import seaborn as sns
plt.figure(figsize=(8, 6))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', center=0,
square=True, linewidths=1, cbar_kws={"shrink": 0.8})
plt.title('Correlation Matrix Heatmap', fontsize=14)
plt.tight_layout()
plt.show()
# Calculate VIF
X = df[['sqft', 'bedrooms', 'bathrooms', 'age']]
X_with_const = sm.add_constant(X)
vif_data = pd.DataFrame()
vif_data["Feature"] = X.columns
vif_data["VIF"] = [variance_inflation_factor(X_with_const.values, i+1)
for i in range(len(X.columns))]
print("\n\nVariance Inflation Factor (VIF):")
print(vif_data)
print("\nInterpretation:")
print("- VIF < 5: Low multicollinearity (acceptable)")
print("- 5 ≤ VIF < 10: Moderate multicollinearity (caution)")
print("- VIF ≥ 10: High multicollinearity (problematic)")
# Remedies for multicollinearity
print("\n\nRemedies for Multicollinearity:")
print("1. Remove highly correlated features")
print("2. Use regularization (Ridge or Lasso regression)")
print("3. Combine correlated features (e.g., feature engineering)")
print("4. Use Principal Component Analysis (PCA)")
# Example: Remove highly correlated feature
X_reduced = df[['sqft', 'age']] # Remove bedrooms and bathrooms
X_reduced_const = sm.add_constant(X_reduced)
vif_reduced = pd.DataFrame()
vif_reduced["Feature"] = X_reduced.columns
vif_reduced["VIF"] = [variance_inflation_factor(X_reduced_const.values, i+1)
for i in range(len(X_reduced.columns))]
print("\n\nVIF after removing correlated features:")
print(vif_reduced)High multicollinearity (VIF ≥ 10) can cause unstable coefficients and inflated standard errors. Use correlation matrices and VIF scores to detect and address multicollinearity issues.
Regularization: Ridge and Lasso Regression
Regularization techniques help prevent overfitting and handle multicollinearity by adding penalty terms to the cost function. Ridge (L2) and Lasso (L1) are common approaches.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge, Lasso, LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split, cross_val_score
# Generate sample data
np.random.seed(42)
n = 100
sqft = np.random.normal(1500, 500, n)
bedrooms = np.random.randint(1, 6, n)
age = np.random.normal(10, 5, n)
price = 100 * sqft + 20000 * bedrooms - 5000 * age + np.random.normal(0, 50000, n)
df = pd.DataFrame({
'sqft': sqft,
'bedrooms': bedrooms,
'age': age,
'price': price
})
X = df[['sqft', 'bedrooms', 'age']]
y = df['price']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Standardize features for regularization
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# 1. Standard Linear Regression
lr = LinearRegression()
lr.fit(X_train_scaled, y_train)
lr_pred = lr.predict(X_test_scaled)
lr_r2 = r2_score(y_test, lr_pred)
lr_mse = mean_squared_error(y_test, lr_pred)
print("Standard Linear Regression:")
print(f"R²: {lr_r2:.4f}")
print(f"MSE: {lr_mse:.2f}")
print(f"Coefficients: {lr.coef_}")
# 2. Ridge Regression (L2 regularization)
# Ridge adds penalty: λ∑βᵢ²
ridge = Ridge(alpha=1.0) # alpha is the regularization strength
ridge.fit(X_train_scaled, y_train)
ridge_pred = ridge.predict(X_test_scaled)
ridge_r2 = r2_score(y_test, ridge_pred)
ridge_mse = mean_squared_error(y_test, ridge_pred)
print("\nRidge Regression (α=1.0):")
print(f"R²: {ridge_r2:.4f}")
print(f"MSE: {ridge_mse:.2f}")
print(f"Coefficients: {ridge.coef_}")
# 3. Lasso Regression (L1 regularization)
# Lasso adds penalty: λ∑|βᵢ|
lasso = Lasso(alpha=1.0)
lasso.fit(X_train_scaled, y_train)
lasso_pred = lasso.predict(X_test_scaled)
lasso_r2 = r2_score(y_test, lasso_pred)
lasso_mse = mean_squared_error(y_test, lasso_pred)
print("\nLasso Regression (α=1.0):")
print(f"R²: {lasso_r2:.4f}")
print(f"MSE: {lasso_mse:.2f}")
print(f"Coefficients: {lasso.coef_}")
# Compare coefficient magnitudes
coef_comparison = pd.DataFrame({
'Feature': X.columns,
'Linear': lr.coef_,
'Ridge': ridge.coef_,
'Lasso': lasso.coef_
})
print("\n\nCoefficient Comparison:")
print(coef_comparison)
# Tune regularization parameter (alpha)
alphas = np.logspace(-4, 2, 100)
ridge_scores = []
lasso_scores = []
for alpha in alphas:
ridge_cv = Ridge(alpha=alpha)
lasso_cv = Lasso(alpha=alpha)
ridge_score = cross_val_score(ridge_cv, X_train_scaled, y_train,
cv=5, scoring='r2').mean()
lasso_score = cross_val_score(lasso_cv, X_train_scaled, y_train,
cv=5, scoring='r2').mean()
ridge_scores.append(ridge_score)
lasso_scores.append(lasso_score)
# Find best alpha
best_ridge_alpha = alphas[np.argmax(ridge_scores)]
best_lasso_alpha = alphas[np.argmax(lasso_scores)]
print(f"\n\nBest Ridge alpha: {best_ridge_alpha:.4f}")
print(f"Best Lasso alpha: {best_lasso_alpha:.4f}")
# Visualize alpha vs R²
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(alphas, ridge_scores, label='Ridge', linewidth=2)
plt.axvline(x=best_ridge_alpha, color='red', linestyle='--',
label=f'Best α={best_ridge_alpha:.4f}')
plt.xscale('log')
plt.xlabel('Alpha (log scale)', fontsize=12)
plt.ylabel('R² Score (CV)', fontsize=12)
plt.title('Ridge Regression: Alpha vs R²', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(1, 2, 2)
plt.plot(alphas, lasso_scores, label='Lasso', linewidth=2, color='orange')
plt.axvline(x=best_lasso_alpha, color='red', linestyle='--',
label=f'Best α={best_lasso_alpha:.4f}')
plt.xscale('log')
plt.xlabel('Alpha (log scale)', fontsize=12)
plt.ylabel('R² Score (CV)', fontsize=12)
plt.title('Lasso Regression: Alpha vs R²', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print("\n\nKey Differences:")
print("- Ridge: Shrinks coefficients toward zero but doesn't eliminate features")
print("- Lasso: Can set coefficients to exactly zero (feature selection)")
print("- Use Ridge when all features are relevant")
print("- Use Lasso when feature selection is desired")Ridge regression (L2) shrinks coefficients but keeps all features, while Lasso regression (L1) can eliminate features by setting coefficients to zero. Elastic Net combines both approaches. Use cross-validation to find optimal regularization strength (α).
Feature Scaling and Transformation
Feature scaling is crucial for regularization and can improve model performance. Log transformations can stabilize variance for skewed targets.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split
# Generate skewed data
np.random.seed(42)
sqft = np.random.lognormal(np.log(1500), 0.5, 100) # Log-normal distribution
price = 100 * sqft + np.random.lognormal(np.log(50000), 0.3, 100)
df = pd.DataFrame({
'sqft': sqft,
'price': price
})
X = df[['sqft']]
y = df['price']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# 1. No scaling
lr_no_scale = LinearRegression()
lr_no_scale.fit(X_train, y_train)
pred_no_scale = lr_no_scale.predict(X_test)
r2_no_scale = r2_score(y_test, pred_no_scale)
# 2. Standard scaling
scaler_standard = StandardScaler()
X_train_std = scaler_standard.fit_transform(X_train)
X_test_std = scaler_standard.transform(X_test)
lr_std = LinearRegression()
lr_std.fit(X_train_std, y_train)
pred_std = lr_std.predict(X_test_std)
r2_std = r2_score(y_test, pred_std)
# 3. Log transformation of target
y_train_log = np.log(y_train)
y_test_log = np.log(y_test)
lr_log = LinearRegression()
lr_log.fit(X_train_std, y_train_log)
pred_log = np.exp(lr_log.predict(X_test_std)) # Transform back
r2_log = r2_score(y_test, pred_log)
print("Model Performance Comparison:")
print(f"No scaling - R²: {r2_no_scale:.4f}")
print(f"Standard scaling - R²: {r2_std:.4f}")
print(f"Log transform + scaling - R²: {r2_log:.4f}")
# Visualize transformations
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
# Original data
axes[0, 0].scatter(X_train, y_train, alpha=0.6, s=50)
axes[0, 0].set_xlabel('Square Feet', fontsize=12)
axes[0, 0].set_ylabel('Price', fontsize=12)
axes[0, 0].set_title('Original Data (Skewed)', fontsize=14)
axes[0, 0].grid(True, alpha=0.3)
# Log-transformed target
axes[0, 1].scatter(X_train, y_train_log, alpha=0.6, s=50)
axes[0, 1].set_xlabel('Square Feet', fontsize=12)
axes[0, 1].set_ylabel('Log(Price)', fontsize=12)
axes[0, 1].set_title('Log-Transformed Target', fontsize=14)
axes[0, 1].grid(True, alpha=0.3)
# Scaled features
axes[1, 0].scatter(X_train_std, y_train, alpha=0.6, s=50)
axes[1, 0].set_xlabel('Square Feet (Standardized)', fontsize=12)
axes[1, 0].set_ylabel('Price', fontsize=12)
axes[1, 0].set_title('Standardized Features', fontsize=14)
axes[1, 0].grid(True, alpha=0.3)
# Both scaled
axes[1, 1].scatter(X_train_std, y_train_log, alpha=0.6, s=50)
axes[1, 1].set_xlabel('Square Feet (Standardized)', fontsize=12)
axes[1, 1].set_ylabel('Log(Price)', fontsize=12)
axes[1, 1].set_title('Both Scaled and Log-Transformed', fontsize=14)
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print("\n\nBest Practices for Transformations:")
print("1. Standardize features when using regularization")
print("2. Log-transform skewed targets to stabilize variance")
print("3. Consider Box-Cox transformation for optimal results")
print("4. Always transform features consistently (fit on train, transform test)")
print("5. Document all transformations for reproducibility")Standardization (z-score) transforms features to have mean 0 and standard deviation 1, making coefficients comparable. Log transformation of skewed targets can stabilize variance and improve model assumptions. Always fit transformers on training data and apply to test data.
Model Evaluation Metrics
Multiple metrics help evaluate linear regression models. Understanding when to use each metric is crucial for proper model assessment.
import numpy as np
from sklearn.metrics import (mean_squared_error, mean_absolute_error,
r2_score, mean_absolute_percentage_error)
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
# Generate sample data
np.random.seed(42)
X = np.random.normal(1500, 500, 100).reshape(-1, 1)
y = 100 * X.flatten() + np.random.normal(0, 50000, 100)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Train model
model = LinearRegression()
model.fit(X_train, y_train)
# Predictions
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)
# Calculate all metrics
metrics_train = {
'MSE': mean_squared_error(y_train, y_train_pred),
'RMSE': np.sqrt(mean_squared_error(y_train, y_train_pred)),
'MAE': mean_absolute_error(y_train, y_train_pred),
'R²': r2_score(y_train, y_train_pred),
}
metrics_test = {
'MSE': mean_squared_error(y_test, y_test_pred),
'RMSE': np.sqrt(mean_squared_error(y_test, y_test_pred)),
'MAE': mean_absolute_error(y_test, y_test_pred),
'R²': r2_score(y_test, y_test_pred),
}
# MAPE (avoid division by zero)
epsilon = 1e-10
mape_train = np.mean(np.abs((y_train - y_train_pred) / (y_train + epsilon))) * 100
mape_test = np.mean(np.abs((y_test - y_test_pred) / (y_test + epsilon))) * 100
print("Training Metrics:")
for metric, value in metrics_train.items():
print(f" {metric}: {value:.4f}")
print(f" MAPE: {mape_train:.2f}%")
print("\nTest Metrics:")
for metric, value in metrics_test.items():
print(f" {metric}: {value:.4f}")
print(f" MAPE: {mape_test:.2f}%")
print("\n\nMetric Interpretations:")
print("- MSE (Mean Squared Error): Penalizes large errors more (in units²)")
print("- RMSE (Root Mean Squared Error): Same units as target, easier to interpret")
print("- MAE (Mean Absolute Error): Average absolute error, robust to outliers")
print("- R² (Coefficient of Determination): Proportion of variance explained (0-1, higher is better)")
print("- MAPE (Mean Absolute Percentage Error): Percentage error, good for comparisons")
print("\n\nWhen to Use Each Metric:")
print("- Use RMSE when large errors are particularly problematic")
print("- Use MAE when all errors are equally important")
print("- Use R² for overall model fit assessment")
print("- Use MAPE for percentage-based business metrics")Different metrics provide different insights. RMSE emphasizes large errors, MAE treats all errors equally, R² measures variance explained, and MAPE provides percentage-based interpretation. Always evaluate on a separate test set to avoid overfitting.
Best Practices
- Split data properly: Always use train/validation/test splits to avoid overfitting.
- Standardize features: Scale features when using regularization or comparing coefficients.
- Check assumptions: Verify linearity, homoscedasticity, normality, and independence through residual analysis.
- Handle multicollinearity: Use VIF scores and correlation matrices to detect and address correlated features.
- Use regularization: Apply Ridge or Lasso when features are numerous or collinear.
- Transform when needed: Log-transform skewed targets to stabilize variance.
- Feature engineering: Create meaningful features through interaction terms, polynomial features, or domain knowledge.
- Cross-validation: Use k-fold cross-validation for robust model evaluation and hyperparameter tuning.
- Interpret coefficients carefully: Consider feature scaling and correlation when interpreting coefficients.
- Visualize results: Use residual plots, Q-Q plots, and prediction plots to diagnose model performance.
- Document assumptions: Clearly state model assumptions and share analysis through
info.studygrid@gmail.com.
Next Steps
Extend linear regression with polynomial features to capture non-linear trends.
Try It Yourself
Exercise 1: Fit a line to hours studied vs score and predict the score for 6 hours.
Show solution
import numpy as np
from sklearn.linear_model import LinearRegression
X = np.array([1, 2, 3, 4, 5]).reshape(-1, 1)
y = np.array([52, 60, 68, 77, 90])
model = LinearRegression().fit(X, y)
print(model.predict([[6]]).round(1)) # ~ 96.5Exercise 2: An R² of 0.2 means what about your model's fit?
Show solution
The model explains only 20% of the variance in the target — a weak fit. The relationship may be non-linear or need more features.
Key Takeaways
- Linear regression fits a straight line minimizing squared error.
- Coefficients show each feature's effect on the target.
- R² measures fit; check residuals and assumptions.
📘 Real-World Deep Dive
Knowing <strong>ML Linear Regression (scikit-learn)</strong> well is what turns scikit-learn from a curiosity into a daily tool — you'll reach for it in nearly every real project.
Real-Life Scenario
An end-to-end usage of ML Linear Regression that you'd actually see in a data pipeline or analytics notebook.
Real-Life Example
import numpy as np
from sklearn.linear_model import LinearRegression
X = np.array([[i] for i in range(1, 11)])
y = np.array([2*i + 1 for i in range(1, 11)], dtype=float)
m = LinearRegression().fit(X, y)
print("slope:", m.coef_[0], "intercept:", m.intercept_)Expected Output
(see source)Common mistakes
fitexpects numeric arrays;OneHotEncoder/LabelEncoderare easy to forget for categorical features.- Calling
predicton a model trained on unscaled data and then scaling inputs at inference time silently degrades accuracy. train_test_split(X, y)requires both arrays;train_test_split(X)for unsupervised learning slips past static checkers.- Treating ML Linear Regression as a black box without reading the docs — the API has subtle defaults that bite when you scale.
🚀 Performance & Best Practices
- Wrap preprocessing + estimator in a
Pipelinesofit/predictstay reproducible. - Use
joblib/picklefor serialising models, not the entire Python state. - For > 100 k rows, switch to
HistGradientBoostingClassifierorcuML. - When working with scikit-learn, prefer vectorised / batched operations over Python loops.
🧪 Try It Yourself
- Reproduce the snippet on a representative slice of your own data.
- Profile the snippet with
cProfileortimeitand find the single biggest improvement. - Generalise the snippet into a small, reusable function you can drop into future projects.