Diagnosing and Evaluating a Multiple Linear Regression Model in R

In this article, we will use the Boston data set to explore some statistical properties using R.

First, we identified and removed the chas, rad, and tax discrete variables by checking the number of unique values in each variable.

Fits a linear regression model predicting medv (median home value):

model <- lm(medv ~ . - chas - rad - tax, data = Boston)

Checking Standard Regression Assumptions

These assumptions underpin the statistical properties of ordinary least squares (OLS) estimation, and violating them can lead to biased estimates, incorrect inferences, or misleading predictions.

• Mean Zero Assumption
  Ensures the model’s errors are not systematically biased. If the mean of residuals is not zero, the model systematically overpredicts or underpredicts the response variable.
• Homoscedasticity (Constant Variance)
  Confirm that residual variance is consistent across all predictor values. Heteroscedasticity (uneven variance) reduces the efficiency of coefficient estimates and invalidates F-tests for model significance
• Normality of Residuals
  Validate that residuals follow a normal distribution. While OLS doesn’t require normality for unbiased estimates, normality ensures valid confidence intervals and p-values in small samples

# ---- CHECKING STANDARD ASSUMPTIONS ----

# 1. Mean Zero Assumption
cat("\n--- Mean Zero Assumption ---\n")
mean_residuals <- mean(residuals(model))
cat("Mean of residuals:", mean_residuals, "\n")

# 2. Homoscedasticity (Constant Variance)
cat("\n--- Homoscedasticity Assumption ---\n")

# Visual check
par(mfrow = c(1, 2))
plot(fitted(model), residuals(model), 
     main = "Residuals vs Fitted Values",
     xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, col = "red", lty = 2)

# Scale-Location plot
plot(fitted(model), sqrt(abs(rstandard(model))), 
     main = "Scale-Location Plot",
     xlab = "Fitted Values", ylab = "√|Standardized Residuals|")
lines(lowess(fitted(model), sqrt(abs(rstandard(model)))), col = "red")

# Formal tests for heteroscedasticity
bp_test <- bptest(model)
cat("Breusch-Pagan test for heteroscedasticity:\n")
print(bp_test)


# 3. Normality of Residuals
cat("\n--- Normality of Residuals Assumption ---\n")

# Visual check
par(mfrow = c(1, 2))
hist(residuals(model), 
     main = "Histogram of Residuals",
     xlab = "Residuals", freq = FALSE)
lines(density(residuals(model)), col = "blue", lwd = 2)
curve(dnorm(x, mean = mean(residuals(model)), sd = sd(residuals(model))), 
      add = TRUE, col = "red", lty = 2, lwd = 2)
legend("topright", c("Density", "Normal"), 
       col = c("blue", "red"), lty = c(1, 2), lwd = 2)

# Q-Q plot
qqnorm(residuals(model))
qqline(residuals(model), col = "red")

# Formal tests for normality
ks_test <- lillie.test(residuals(model))

cat("Kolmogorov-Smirnov test for normality:\n")
print(ks_test)

Outlier Detection

• Outliers in Predictors (Leverage)

Leverage measures how far an observation’s predictor values deviate from the average predictor values in a dataset. It identifies observations with unusual predictor values that may disproportionately influence regression results.
Mathematically, it is the diagonal element hii of the hat matrix H, defined as:

where X is the design matrix of predictor variables. High-leverage points (hii>2p/n, where p = number of predictors, n = sample size) can distort coefficient estimates and predictions.

# ---- CHECKING FOR OUTLIERS IN PREDICTORS ----
cat("\n--- Outliers in Predictors ---\n")

# Calculate leverage (hat values)
leverage <- hatvalues(model)
high_leverage <- which(leverage > 2 * mean(leverage))
cat("Observations with high leverage (> 2 * mean):\n")
print(high_leverage)

# Plot leverage
par(mfrow = c(1, 1))
plot(leverage, 
     main = "Leverage Points",
     ylab = "Leverage", xlab = "Observation Number")
abline(h = 2 * mean(leverage), col = "red", lty = 2)
text(high_leverage, leverage[high_leverage], labels = high_leverage, pos = 3)

# Examine the most significant outlier in predictors
if (length(high_leverage) > 0) {
  most_extreme <- high_leverage[which.max(leverage[high_leverage])]
  cat("\nMost significant outlier in predictors (observation", most_extreme, "):\n")
  print(Boston[most_extreme, ])
  
  # Compare with average values
  cat("\nCompared to average values:\n")
  comparison <- rbind(Boston[most_extreme, ], colMeans(Boston))
  rownames(comparison) <- c("Outlier", "Average")
  print(comparison)
}

• Outliers in Response (Standardized Residuals)

Outliers in response are data points where the observed value of the dependent variable (Y) deviates significantly from the value predicted by the regression model. Detect observations where the model fails to explain the response variable adequately.

Get Mia Lai’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

Standardized Residuals are calculated as:

# ---- CHECKING FOR OUTLIERS IN RESPONSE ----
cat("\n--- Outliers in Response ---\n")

# Calculate standardized residuals
std_resid <- rstandard(model)
outliers_response <- which(abs(std_resid) > 2)
cat("Observations with standardized residuals > 2:\n")
print(outliers_response)

# Plot standardized residuals
par(mfrow = c(1, 1))
plot(std_resid, 
     main = "Standardized Residuals",
     ylab = "Standardized Residuals", xlab = "Observation Number")
abline(h = c(-2, 2), col = "red", lty = 2)
text(outliers_response, std_resid[outliers_response], labels = outliers_response, pos = 3)

# Examine the most significant outlier in response
if (length(outliers_response) > 0) {
  most_extreme_response <- outliers_response[which.max(abs(std_resid[outliers_response]))]
  cat("\nMost significant outlier in response (observation", most_extreme_response, "):\n")
  print(Boston[most_extreme_response, ])
  cat("Standardized residual:", std_resid[most_extreme_response], "\n")
  
  # Compare with predicted value
  cat("Actual medv:", Boston[most_extreme_response, "medv"], "\n")
  cat("Predicted medv:", fitted(model)[most_extreme_response], "\n")
  cat("Difference:", Boston[most_extreme_response, "medv"] - fitted(model)[most_extreme_response], "\n")
}

Influential Observations

• Calculates Cook’s distance, DFFITS, and DFBETAS to identify points that strongly influence the model fit.
• Flags observations exceeding typical thresholds (e.g., Cook’s distance > 4/n).
• Plots Cook’s distance and highlights influential points.
• Compares model coefficients with and without the most influential observation to assess its impact.

# ---- CHECKING FOR INFLUENTIAL OBSERVATIONS ----
cat("\n--- Influential Observations ---\n")

# Calculate Cook's distance
cooks_d <- cooks.distance(model)
influential <- which(cooks_d > 4/nrow(Boston))
cat("Influential observations (Cook's distance > 4/n):\n")
print(influential)

# Plot Cook's distance
par(mfrow = c(1, 1))
plot(cooks_d, 
     main = "Cook's Distance",
     ylab = "Cook's Distance", xlab = "Observation Number")
abline(h = 4/nrow(Boston), col = "red", lty = 2)
text(influential, cooks_d[influential], labels = influential, pos = 3)

# DFFITS
dffits_val <- dffits(model)
dffits_threshold <- 2 * sqrt(length(coef(model))/nrow(Boston))
influential_dffits <- which(abs(dffits_val) > dffits_threshold)
cat("\nInfluential observations (DFFITS > threshold):\n")
print(influential_dffits)

# DFBETAS
dfbetas_val <- dfbetas(model)
dfbetas_threshold <- 2/sqrt(nrow(Boston))
influential_dfbetas <- which(apply(abs(dfbetas_val), 1, max) > dfbetas_threshold)
cat("\nInfluential observations (DFBETAS > threshold):\n")
print(influential_dfbetas)

# Examine the most influential observation
if (length(influential) > 0) {
  most_influential <- influential[which.max(cooks_d[influential])]
  cat("\nMost influential observation (", most_influential, "):\n")
  print(Boston[most_influential, ])
  cat("Cook's distance:", cooks_d[most_influential], "\n")
  
  # Compare model with and without the most influential observation
  model_without <- lm(medv ~ . - chas - rad, data = Boston[-most_influential, ])
  cat("\nCoefficients comparison (with vs. without influential point):\n")
  coef_comparison <- cbind(coef(model), coef(model_without), 
                           (coef(model) - coef(model_without))/coef(model) * 100)
  colnames(coef_comparison) <- c("Original", "Without Influential", "% Change")
  print(coef_comparison)
}

Intersection Method

• Identifies points that are both high-leverage and have large residuals (potentially most problematic).
• Plots leverage vs. absolute standardized residuals, highlighting these influential points.
• Compares model coefficients with and without these points.

# ---- MODIFIED INFLUENTIAL OBSERVATION CRITERIA ----
cat("\n--- Influential Observations (Intersection Method) ---\n")

# 1. Identify predictor outliers (high leverage)
leverage <- hatvalues(model)
high_leverage <- leverage > 2 * mean(leverage)

# 2. Identify response outliers (large residuals)
std_resid <- rstandard(model)
large_residuals <- abs(std_resid) > 2

# 3. Find intersection of both conditions
influential_intersection <- which(high_leverage & large_residuals)
cat("Influential observations (High Leverage AND Large Residuals):\n")
print(influential_intersection)

# 4. Visualize the relationship
par(mfrow = c(1, 1))
plot(leverage, abs(std_resid),
     main = "Leverage vs. Absolute Standardized Residuals",
     xlab = "Leverage", ylab = "|Standardized Residuals|",
     col = ifelse(high_leverage & large_residuals, "red", "black"))
abline(v = 2*mean(leverage), col = "blue", lty = 2)
abline(h = 2, col = "blue", lty = 2)
legend("topright", 
       legend = c("Influential", "Thresholds"),
       col = c("red", "blue"), pch = c(1, NA), lty = c(NA, 2))

# 5. Examine characteristics of influential points
if (length(influential_intersection) > 0) {
  cat("\nCharacteristics of influential observations:\n")
  print(Boston[influential_intersection, ])
  
  # Compare with original model
  cat("\nOriginal model coefficients:\n")
  print(coef(model))
  
  # Model without influential points
  model_clean <- update(model, subset = -influential_intersection)
  cat("\nModel coefficients without influential points:\n")
  print(coef(model_clean))
}

Multicollinearity Diagnostics

• Computes the correlation matrix among predictors and highlights pairs with high correlation (> 0.7).
• Calculates Variance Inflation Factors (VIF) for each predictor; VIF > 5 suggests problematic multicollinearity, > 10 severe.
• Calculates the condition number of the correlation matrix; values > 30 indicate multicollinearity.
• Plots VIFs with reference lines for interpretation.

# ---- CHECKING OR MULTICOLLINEARITY ----
cat("\n--- Multicollinearity ---\n")

# Correlation matrix
predictors <- Boston[, !(names(Boston) %in% c("medv", "chas", "rad"))]
cor_matrix <- cor(predictors)
cat("Correlation matrix:\n")
print(round(cor_matrix, 2))

# Identify high correlations
high_cors <- which(abs(cor_matrix) > 0.7 & abs(cor_matrix) < 1, arr.ind = TRUE)
if (nrow(high_cors) > 0) {
  cat("\nHigh correlations (> 0.7):\n")
  for (i in 1:nrow(high_cors)) {
    if (high_cors[i, 1] < high_cors[i, 2]) {  # Avoid duplicates
      cat(colnames(predictors)[high_cors[i, 1]], "and", 
          colnames(predictors)[high_cors[i, 2]], ":", 
          cor_matrix[high_cors[i, 1], high_cors[i, 2]], "\n")
    }
  }
}

# Variance Inflation Factors (VIF)
vif_values <- vif(model)
cat("\nVariance Inflation Factors:\n")
print(vif_values)
cat("\nVIF > 5 indicates problematic multicollinearity\n")
cat("VIF > 10 indicates severe multicollinearity\n")

# Condition number
eigen_values <- eigen(cor_matrix)$values
condition_number <- sqrt(max(eigen_values) / min(eigen_values))
cat("\nCondition Number:", condition_number, "\n")
cat("Condition Number > 30 indicates multicollinearity\n")

# Plot VIF values
par(mfrow = c(1, 1))
barplot(vif_values, main = "Variance Inflation Factors", 
        horiz = TRUE, col = "steelblue", las = 1)
abline(v = c(5, 10), col = c("orange", "red"), lty = 2)
text(x = c(5, 10), y = par("usr")[3] - 1, labels = c("Problematic", "Severe"), 
     col = c("orange", "red"), pos = 4)