# Foundations R Regression Correlation And Diagnostics

> Fit/diagnose R lm() regression and compute Pearson/Spearman/Kendall correlation with Fisher z CIs. Use for lm() regression, residual diagnostics, correlation CIs, or t.test/shapiro.test checks.

- Skill: `pavel-kravchenko/foundations-r-regression-correlation-and-diagnostics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/foundations-r-regression-correlation-and-diagnostics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/foundations-r-regression-correlation-and-diagnostics/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

| Situation | Function | Distribution |
|---|---|---|
| Mean, σ unknown (usual case) | `qt(1-α/2, df=n-1)` | Student-t |
| Variance σ² | `qchisq(...)` | Chi-squared |
| Count / Poisson mean, large n | `qnorm(1-α/2)` | Normal asymptotic |
| Proportion, large n | `qnorm(1-α/2)` | Normal asymptotic |

```r
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.

```r
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

```r
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.

```r
#' 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)

```r
#' 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².

```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

```r
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.