Unit 5 — Regression Modeling in R

Topics Covered

1. Correlation Analysis

1.1 Karl Pearson Correlation Coefficient

Measures the strength of linear association between two numeric variables; range −1 to +1.

# Sample data
hours <- c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11)
marks <- c(50, 55, 62, 66, 72, 75, 80, 85, 90, 95)

cor(hours, marks)                          # Pearson (default)
# [1] 0.9985

# Same explicitly
cor(hours, marks, method = "pearson")

1.2 Spearman Rank Correlation

Non-parametric; based on ranks; captures any monotonic relationship (linear or not). Robust to outliers.

cor(hours, marks, method = "spearman")
# [1] 1.000  (perfect monotone increase here)

1.3 Significance Test — cor.test()

cor.test(hours, marks)
# Returns: t-statistic, df, p-value, 95 % CI, sample estimate

cor.test(hours, marks, method = "spearman")

1.4 Correlation Matrix for Multiple Variables

data(iris)
round(cor(iris[, 1:4]), 3)
#              Sepal.Length Sepal.Width Petal.Length Petal.Width
# Sepal.Length        1.000      -0.118        0.872       0.818
# Sepal.Width        -0.118       1.000       -0.428      -0.366
# Petal.Length        0.872      -0.428        1.000       0.963
# Petal.Width         0.818      -0.366        0.963       1.000

# Visualise it
# install.packages("corrplot")
library(corrplot)
corrplot(cor(iris[, 1:4]), method = "color",
         addCoef.col = "black", tl.col = "black")

1.5 Karl Pearson Coefficient from Scratch

n <- length(hours)
num <- n * sum(hours * marks) - sum(hours) * sum(marks)
den <- sqrt((n * sum(hours^2) - sum(hours)^2) *
            (n * sum(marks^2) - sum(marks)^2))
num / den                                  # equals cor(hours, marks)

1.6 Spearman from Scratch (using ranks)

rx <- rank(hours)
ry <- rank(marks)
cor(rx, ry)                                # equals Spearman ρ
1 - 6 * sum((rx - ry)^2) / (n * (n^2 - 1))  # classical formula

cor.test(iris$Sepal.Length, iris$Petal.Length)
# t = 21.65, p < 2.2e-16,  r = 0.872, 95 % CI (0.827, 0.906)

x <- c(80, 85, 85, 70, 90)        # tie at 85
y <- c(70, 80, 90, 60, 95)
cor(x, y, method = "spearman")    # 0.975
cor.test(x, y, method = "spearman")

2. Simple Linear Regression with lm()

2.1 Fitting the Model

fit <- lm(marks ~ hours)             # formula: Y ~ X
fit                                  # short print: coefficients only

Fitted line: \(\hat Y = 41.09 + 4.91\, X\).

2.2 Full Summary

summary(fit)

• Estimate — coefficient values.
• Std. Error — standard error of each coefficient.
• t value — Estimate / Std. Error.
• Pr(>|t|) — p-value testing H₀: coefficient = 0.
• R-squared — proportion of variation explained.
• F-statistic — overall model significance.

2.3 Extracting Components

coef(fit)             # coefficients
fitted(fit)           # predicted Y for the training data
resid(fit)            # residuals
fit$model             # data used
confint(fit)          # 95 % CIs for coefficients
anova(fit)            # ANOVA table

# Specific bits
summary(fit)$r.squared     # 0.997
summary(fit)$sigma         # residual SE
deviance(fit)              # SSE (sum of squared residuals)
AIC(fit); BIC(fit)         # information criteria

2.4 Plotting the Regression Line

plot(hours, marks, pch = 19, col = "navy",
     main = "Marks vs Study Hours",
     xlab = "Hours", ylab = "Marks")
abline(fit, col = "red", lwd = 2)
text(3, 90, paste0("R² = ", round(summary(fit)$r.squared, 3)),
     pos = 4, col = "darkred")

2.5 Prediction

# Predict for new X
newx <- data.frame(hours = c(5, 7.5, 12))
predict(fit, newx)                           # point estimates
predict(fit, newx, interval = "confidence")  # CI for mean response
predict(fit, newx, interval = "prediction")  # CI for individual response

Note the prediction intervals are always wider than the confidence intervals — an individual response varies more than a mean.

Always be cautious extrapolating outside the observed range of X.

3. Multiple Linear Regression (preview)

data(mtcars)
fit2 <- lm(mpg ~ wt + hp + cyl, data = mtcars)
summary(fit2)

Just add predictors with +. Use * for interactions; poly(x, 2) for polynomial terms.

lm(mpg ~ wt + I(wt^2), data = mtcars)        # quadratic in wt
lm(mpg ~ wt * hp, data = mtcars)             # main effects + wt:hp interaction

4. Residual Diagnostics

par(mfrow = c(2, 2))
plot(fit)
par(mfrow = c(1, 1))

Four standard panels (recall Unit 3):

1. Residuals vs Fitted — linearity / homoscedasticity.
2. Normal Q-Q — normality of residuals.
3. Scale–Location — checks equal variance.
4. Residuals vs Leverage — influential observations (Cook's distance).

Formal Checks

shapiro.test(resid(fit))             # Normality
# install.packages("car")
library(car)
ncvTest(fit)                         # Heteroscedasticity
durbinWatsonTest(fit)                # Auto-correlation of residuals
vif(fit2)                            # Multicollinearity (multiple regression)

Cook's Distance — Influential Points

cooks <- cooks.distance(fit2)
plot(cooks, type = "h", main = "Cook's Distance", ylab = "")
abline(h = 4 / nrow(mtcars), col = "red", lty = 2)
# Rule of thumb: distance > 4/n is suspect

5. Comparing Pearson vs Spearman Visually

# Generate non-linear monotone data
x <- 1:50
y <- exp(x/15) + rnorm(50, 0, 5)

cor(x, y)                                 # Pearson — weaker (linear only)
cor(x, y, method = "spearman")            # Spearman — closer to 1

plot(x, y, pch = 19, col = "navy",
     main = "Pearson sees the noise, Spearman sees the order")
abline(lm(y ~ x), col = "red", lwd = 2)   # poor linear fit
lines(x, exp(x/15), col = "darkgreen", lwd = 2, lty = 2)  # true curve

6. Standard Regression Workflow

1. Explore the data: summary(), cor(), scatter plots.
2. Fit the model with lm().
3. Examine summary() — coefficient p-values, R², overall F.
4. Diagnose residuals with plot(fit).
5. Refine — drop non-significant predictors, transform variables (log(), sqrt()), add interactions or polynomial terms.
6. Validate — predict on held-out data; compute RMSE / MAE.
7. Report — equation, R², residual SE, key p-values, plot the line + CI.

# Train-test split & RMSE
set.seed(1)
idx   <- sample(seq_len(nrow(mtcars)), 0.7 * nrow(mtcars))
train <- mtcars[idx, ]
test  <- mtcars[-idx, ]

fit <- lm(mpg ~ wt + hp, data = train)
pred <- predict(fit, newdata = test)
rmse <- sqrt(mean((test$mpg - pred)^2))
rmse

fit <- lm(mpg ~ wt, data = mtcars)
summary(fit)
plot(mtcars$wt, mtcars$mpg, pch = 19, col = "darkblue",
     main = "MPG vs Weight (mtcars)")
abline(fit, col = "red", lwd = 2)

new <- data.frame(wt = c(2.5, 3.5))
predict(fit, new, interval = "prediction")  # predict MPG for new weights

x <- 1:50
y <- 0.5 * x^2 + 3 * x + rnorm(50, 0, 20)

fit_lin  <- lm(y ~ x)
fit_quad <- lm(y ~ poly(x, 2))
summary(fit_lin)$adj.r.squared       # lower
summary(fit_quad)$adj.r.squared      # higher

plot(x, y, pch = 19, col = "navy")
abline(fit_lin, col = "red", lwd = 2)
lines(x, predict(fit_quad), col = "green", lwd = 2)
legend("topleft", c("Linear","Quadratic"),
       col = c("red","green"), lty = 1, lwd = 2, bty = "n")

7. Useful Helper Functions

Key Take-aways

• cor() for correlation; default Pearson; use method = "spearman" for rank.
• cor.test() tests H₀: ρ = 0 and gives a CI.
• lm(y ~ x) fits a linear model; summary() reports coefficients, R², F-test.
• predict() with interval = "confidence" or "prediction".
• plot(fit) gives four standard residual diagnostics.
• Spearman is robust to non-linearity and outliers; use when linearity / normality fail.