Foundations R Regression Correlation And Diagnostics

R Regression, Correlation, and Diagnostics

When to Use

• Fitting lm() simple/multiple/polynomial regression and need to check whether the fit is trustworthy (not just R²).
• Choosing between Pearson, Spearman, or Kendall correlation for a pair of continuous variables.
• Building a confidence interval for a mean, variance, proportion, or count-rate, or reporting one alongside a t.test/cor.test.
• Deciding between a parametric test (t-test, Pearson) and its nonparametric counterpart (Wilcoxon, Spearman) based on a normality check.
• Running one-way ANOVA via lm + anova() and validating it against Kruskal-Wallis on the same groups.

Version Compatibility

Base R ≥ 4.0 (stats package: lm, cor.test, t.test, shapiro.test, qqnorm are all base/stats, no install needed). Optional: nortest ≥ 1.0 (pearson.test), exactRankTests ≥ 0.8 (wilcox.exact, exact p-values for small n with ties). Notebook source: Course/Tier_0_Computational_Foundations/08_Advanced_R_Statistics/02_r_regression_correlation_and_diagnostics.ipynb.

Prerequisites

• install.packages(c("nortest", "exactRankTests")) if using those two functions.
• Comfortable with R data frames (read.table, $, [[) and basic hypothesis-testing concepts (H0/H1, p-value, α).
• Related skill: foundations-r-hypothesis-testing-and-nonparametrics for the nonparametric tests referenced below.

Key Patterns

CI selection

n <- length(x); xbar <- mean(x); s <- sd(x); alpha <- 0.05
q <- qt(1 - alpha / 2, df = n - 1)
ci <- c(xbar - q * s / sqrt(n), xbar + q * s / sqrt(n))
t.test(x, conf.level = 1 - alpha)   # equivalent, also gives the CI directly

Two-sample comparison workflow

Always test variance equality before picking var.equal in t.test; fall back to Wilcoxon if normality is doubtful.

var.test(x, y, ratio = 1, alternative = "two.sided")            # F-test on variances first
t.test(x, y, alternative = "greater", paired = FALSE,
       var.equal = FALSE)                                        # Welch's t-test (safe default)

library(exactRankTests)
wilcox.exact(x, y, paired = FALSE, alternative = "greater", exact = TRUE)  # nonparametric check

Normality testing

shapiro.test(x)                       # most powerful for n <= 5000
library(nortest)
pearson.test(x, adjust = FALSE)       # chi-sq goodness-of-fit; true p-value between adjust=F/T
qqnorm(x); qqline(x, col = "red")     # always pair a numeric test with the visual check

Goal: decide, before running cor.test, whether Pearson is valid or Spearman/Kendall is safer. Approach: Shapiro-Wilk both variables; if either rejects normality or the scatter looks nonlinear, use rank-based correlation.

#' Pick a correlation method from a quick normality + scatter check
#' @param x,y numeric vectors of equal length
#' @return the fitted htest object from cor.test
choose_correlation <- function(x, y, alpha = 0.05) {
  p_x <- shapiro.test(x)$p.value
  p_y <- shapiro.test(y)$p.value
  method <- if (p_x > alpha && p_y > alpha) "pearson" else "spearman"
  cor.test(x, y, method = method, conf.level = 1 - alpha)
}

Fisher z-transform CI for Pearson r (manual)

#' 95% CI for Pearson's r via the Fisher z-transform
#' @param x,y numeric vectors; alpha significance level
fisher_r_ci <- function(x, y, alpha = 0.05) {
  n  <- length(x)
  r  <- cor(x, y)
  z1 <- 0.5 * log((1 + r) / (1 - r))          # Fisher z of the observed r
  q  <- qnorm(1 - alpha / 2)
  z_lo <- z1 - q / sqrt(n - 3)
  z_hi <- z1 + q / sqrt(n - 3)
  r_ci <- (exp(2 * c(z_lo, z_hi)) - 1) / (exp(2 * c(z_lo, z_hi)) + 1)  # back-transform
  list(r = r, ci = r_ci)
}

Regression fit + residual diagnostics

Goal: fit medv ~ lstat on the Boston housing data and verify the linear model is adequate, not just report R². Approach: fit with lm, inspect summary/confint, then check residuals — a curve or funnel in Residuals-vs-Fitted means the model is misspecified even with high R².

boston_lm <- lm(medv ~ lstat, data = Boston_data)
summary(boston_lm)          # coefficients, R^2, F-statistic
confint(boston_lm)          # 95% CI for slope/intercept

par(mfrow = c(2, 2))
plot(boston_lm)             # Residuals vs Fitted, Q-Q, Scale-Location, Leverage
par(mfrow = c(1, 1))

shapiro.test(residuals(boston_lm))   # residual normality

# Nonlinearity? add a quadratic or orthogonal-polynomial term
lstat2 <- Boston_data$lstat^2
boston_lm_2 <- lm(medv ~ lstat + lstat2, data = Boston_data)
boston_lm_poly7 <- lm(medv ~ poly(lstat, 7), data = Boston_data)  # higher-degree fit
anova(boston_lm, boston_lm_2)        # nested-model F-test: does the quadratic term help?

ANOVA via lm vs. Kruskal-Wallis

lm_result <- lm(revenue ~ store_id, data = pharmacy_data)
anova(lm_result)                                  # parametric one-way ANOVA
kruskal.test(pharmacy_data$revenue, pharmacy_data$store_id)  # nonparametric cross-check
shapiro.test(residuals(lm_result))                # verify ANOVA's normality assumption

Pitfalls

• R² alone is misleading: high R² does not mean the model is correct or causal — always check plot(model), especially Residuals vs. Fitted (pattern = misspecification; random scatter around zero = good).
• Pearson vs. Spearman: Pearson measures linear association and assumes approximately normal data; Spearman/Kendall measure monotonic association via ranks — use for skewed data (counts, survival times) or when linearity is uncertain.
• Confounders in multiple regression: adding a variable changes every other coefficient — each coefficient means "effect holding all others constant," not a marginal effect.
• var.test before t.test: always check variance equality first; var.equal = TRUE on unequal variances inflates Type I error.
• Bonferroni is conservative: for many correlated tests (genomic screens), prefer Benjamini-Hochberg (p.adjust(method = "BH")) over Bonferroni.
• Polynomial overfitting: a degree-7 poly() fit can track noise, not signal — compare nested models with anova(), don't just chase R².

See Also

• foundations-r-hypothesis-testing-and-nonparametrics — Wilcoxon, sign test, Kruskal-Wallis with Dunn post-hoc.
• bio-experimental-design-multiple-testing — FDR/BH correction for many simultaneous tests.
• statistical-analysis — general Python-side statistical workflow.
• bio-differential-expression-deseq2-basics — regression-style modeling (GLMs) for count-based omics data.