Assumption Checklist Template (Statistics, Diagnostics, R and Python)

Assumption Checklist: [Your study question or course]

A working checklist for the assumptions behind the methods you will use most: t-tests and ANOVA, linear regression, generalised linear models (logistic and Poisson), chi-square tests, correlation, and repeated-measures and mixed models. For each method family it lists the assumption, the diagnostic plot and test that checks it, the R and Python code (run on one shared data set, with the real output), how to judge the result, and what to do if the check fails. It is for students learning why assumptions matter, master's and PhD students who must defend their analysis to a supervisor or examiner, and analysts who want a repeatable checking routine. It speeds up learning by using one small data set for every family, so you see the same outlier and the same design choices from several angles, and by teaching a judgement rule (plot first, size of the problem second, test p-value last). The worked examples use fictional data; replace them with your own. Source basis: the NIST/SEMATECH e-Handbook of Statistical Methods, Box and Cox (1964), Cook (1977), Nelder and Wedderburn (1972), and the R and SciPy/statsmodels documentation, written up in original words and examples.

Document control

Why assumptions matter and how to judge them

Every statistical test is a calculation that gives the right answer if the data behave in a certain way. The p-value of a t-test, for example, is computed as if the sample means followed a t distribution; that is approximately true when the observations are independent and the means are close to normal. When an assumption fails, the method does not crash: it quietly gives a p-value or confidence interval that is too small, too large, or about the wrong quantity. Checking assumptions is how you find out whether the numbers on your screen mean what you think they mean.

Not all assumptions are equal. A useful order of importance, from most to least serious when violated:

1. Independence and correct design. Clustered or repeated observations analysed as independent give standard errors that are too small. No amount of sample size fixes this, and no plot of the outcome reveals it: you check it from how the data were collected.
2. The right model form. Linearity in regression, the right link function in a GLM, the right variance function (Poisson vs negative binomial). A wrong form biases the estimates themselves.
3. Influential points and outliers. One unusual observation can drive a whole result.
4. Equal variances. Matters most when group sizes are unequal; the Welch versions of the t-test and ANOVA remove most of the problem.
5. Normality of residuals. Matters least with moderate or large samples, because the central limit theorem makes the sampling distribution of means and coefficients close to normal. It matters for prediction intervals and with very small samples.

The judgement rule: plot first, size second, test last. Look at a plot that shows the assumption directly (a QQ plot, residuals against fitted values, a dot plot per group). Ask how big the departure is and whether it could change the conclusion. Use a formal test (Shapiro-Wilk, Levene, Breusch-Pagan) as a supporting number, never as an automatic switch.

Why formal normality tests mislead

A normality test answers "is there evidence that this distribution is not exactly normal?". That is the wrong question. No real variable is exactly normal; the useful question is "is the departure large enough to distort my analysis?". Formal tests answer the wrong question badly in both directions, because their power depends on n:

• With small samples, where non-normality matters most, the test has little power and often misses clear skew.
• With large samples, where the central limit theorem protects the t-test and regression, the test detects tiny, harmless departures and returns a very small p-value.

The simulation below shows both effects. It draws 1,000 samples of n = 15 from a clearly right-skewed (lognormal) distribution, and 200 samples of n = 5,000 from a t distribution with 10 df, which looks almost exactly normal.

set.seed(59)
small <- replicate(1000, shapiro.test(rlnorm(15, 0, 0.6))$p.value < 0.05)
big   <- replicate(200,  shapiro.test(rt(5000, df = 10))$p.value < 0.05)
cat("n = 15, skewed (lognormal): rejected", mean(small), "\n")
cat("n = 5000, near-normal (t, df = 10): rejected", mean(big), "\n")

n = 15, skewed (lognormal): rejected 0.508
n = 5000, near-normal (t, df = 10): rejected 1

import numpy as np
from scipy import stats

rng = np.random.default_rng(59)
small = np.mean([stats.shapiro(rng.lognormal(0, 0.6, 15)).pvalue < 0.05 for _ in range(1000)])
big = np.mean([stats.shapiro(rng.standard_t(10, 5000)).pvalue < 0.05 for _ in range(200)])
print("n = 15, skewed (lognormal): rejected", small)
print("n = 5000, near-normal (t, df = 10): rejected", big)

n = 15, skewed (lognormal): rejected 0.492
n = 5000, near-normal (t, df = 10): rejected 1.0

The R and Python numbers differ slightly because the two languages use different random number generators; the pattern is the same. At n = 15 the test missed obvious skew about half the time. At n = 5,000 it rejected a nearly normal distribution every single time. A rule like "if Shapiro-Wilk p < .05, switch to a rank test" therefore switches for the wrong reasons in large samples and fails to switch in small ones.

The shared data set

Twenty-four Lakemoor University (fictional) students in an introductory statistics course were taught in one of three formats (lecture, flipped classroom, online), eight per format. For each student we have weekly study hours, the final exam score (0 to 100), whether they passed (score of 60 or more), the number of library visits in the term (a count), and three quiz scores (out of 10) in weeks 2, 6 and 10. Student 24 is unusual: an online student with modest study hours and a very high exam score. Keep an eye on student 24 in every section.

# Lakemoor University study-habits data (fictional), 24 students
# install.packages(c("car", "lme4", "MASS"))   # once, if not installed
d <- data.frame(
  id     = 1:24,
  format = factor(rep(c("lecture", "flipped", "online"), each = 8),
                  levels = c("lecture", "flipped", "online")),
  hours  = c(4,6,5,8,3,7,5,6,  7,9,6,10,8,5,9,7,  3,5,4,6,2,7,4,5),
  score  = c(56,55,57,66,48,57,62,61,  62,82,68,76,75,57,69,71,  50,51,58,54,42,58,53,78),
  visits = c(0,9,1,15,0,2,6,4,  3,18,1,25,12,0,8,10,  0,7,0,11,0,15,1,2),
  q1 = c(5,6,5,7,4,6,5,6,  6,7,5,8,6,5,7,6,  4,5,5,6,3,6,4,5),
  q2 = c(6,6,7,8,5,6,6,7,  6,8,7,9,8,6,7,7,  5,4,6,7,4,5,5,6),
  q3 = c(7,8,6,8,6,8,7,7,  8,9,7,9,8,7,9,8,  5,6,6,6,5,7,5,7)
)
d$pass <- as.integer(d$score >= 60)

import numpy as np, pandas as pd
from scipy import stats
import statsmodels.api as sm
import statsmodels.formula.api as smf

d = pd.DataFrame({
    "id": range(1, 25),
    "format": ["lecture"] * 8 + ["flipped"] * 8 + ["online"] * 8,
    "hours": [4,6,5,8,3,7,5,6, 7,9,6,10,8,5,9,7, 3,5,4,6,2,7,4,5],
    "score": [56,55,57,66,48,57,62,61, 62,82,68,76,75,57,69,71, 50,51,58,54,42,58,53,78],
    "visits": [0,9,1,15,0,2,6,4, 3,18,1,25,12,0,8,10, 0,7,0,11,0,15,1,2],
    "q1": [5,6,5,7,4,6,5,6, 6,7,5,8,6,5,7,6, 4,5,5,6,3,6,4,5],
    "q2": [6,6,7,8,5,6,6,7, 6,8,7,9,8,6,7,7, 5,4,6,7,4,5,5,6],
    "q3": [7,8,6,8,6,8,7,7, 8,9,7,9,8,7,9,8, 5,6,6,6,5,7,5,7],
})
d["format"] = pd.Categorical(d["format"], ["lecture", "flipped", "online"])
d["pass_"] = (d["score"] >= 60).astype(int)   # "pass" is a Python keyword

Outputs in this template are trimmed for space: R test headers are condensed to one line and blank lines are removed, but every number is copied unchanged from a real run with the software versions in Document control. Blocks that were shortened start with "(output abridged)". Later code blocks assume d (and, in Python, the imports above) already exist.

Overview: assumptions by method family

Family 1: t-tests and ANOVA

The assumptions. Observations are independent within and between groups; there are no extreme outliers; residuals are approximately normal (less important as group sizes grow); and, for the classic F test and the pooled t-test only, the groups have equal variances. Welch's versions drop the equal-variance assumption.

The diagnostics. A dot plot or box plot per group shows outliers and spread at once. A QQ plot of the residuals shows normality. The ratio of the largest to the smallest group SD is a quick spread check: under about 2 the classic F test is usually fine with equal group sizes. Levene's test (median-centred, also called Brown-Forsythe) is the supporting test for equal variances; it is far less sensitive to non-normality than Bartlett's test.

fit_aov <- aov(score ~ format, data = d)
summary(fit_aov)
tapply(d$score, d$format, sd)
shapiro.test(residuals(fit_aov))
car::leveneTest(score ~ format, data = d)   # median-centred by default
oneway.test(score ~ format, data = d)       # Welch ANOVA
kruskal.test(score ~ format, data = d)
# plot(fit_aov, which = 2)                  # QQ plot of residuals

(output abridged)
            Df Sum Sq Mean Sq F value  Pr(>F)   
format       2  974.3   487.2   7.269 0.00399 **
Residuals   21 1407.5    67.0                   
  lecture   flipped    online 
 5.391793  7.964206 10.419761 
Shapiro-Wilk normality test: W = 0.94704, p-value = 0.2336
Levene's Test (center = median): F(2, 21) = 0.511, Pr(>F) = 0.6072
Welch one-way: F = 7.2758, num df = 2.000, denom df = 13.052, p-value = 0.007537
Kruskal-Wallis chi-squared = 9.5537, df = 2, p-value = 0.008423

from statsmodels.stats.oneway import anova_oneway

groups = [g["score"].to_numpy() for _, g in d.groupby("format", observed=True)]
fit_aov = smf.ols("score ~ format", data=d).fit()
print(sm.stats.anova_lm(fit_aov, typ=2).round(4))
print(stats.shapiro(fit_aov.resid))
print(stats.levene(*groups, center="median"))
print(anova_oneway(groups, use_var="unequal"))
print(stats.kruskal(*groups))

(output abridged)
             sum_sq    df       F  PR(>F)
format     974.3333   2.0  7.2686   0.004
Residual  1407.5000  21.0     NaN     NaN
ShapiroResult(statistic=np.float64(0.9470422320797904), pvalue=np.float64(0.23357208958737835))
LeveneResult(statistic=np.float64(0.5109851169383416), pvalue=np.float64(0.6071734827056823))
AnovaResult(statistic=np.float64(7.275783172448055), pvalue=np.float64(0.0075367401472465085), df=(2.0, np.float64(13.051711791153103)), ...)
KruskalResult(statistic=np.float64(9.553672624237143), pvalue=np.float64(0.008422603349467747))

How to judge. The group SDs are 5.4, 8.0 and 10.4, a ratio of 1.93: borderline. Levene's test does not flag it (p = .607), but with eight students per group it has little power, so a non-significant result is not reassurance. The online group's large SD is driven by student 24 (score 78 against a group mean of 55.5). The residuals look acceptable as a whole (Shapiro-Wilk p = .234), but a dot plot would show student 24 sitting far above the rest of the online group. The classic F (F(2, 21) = 7.27, p = .004), Welch's F (F(2, 13.05) = 7.28, p = .008) and Kruskal-Wallis (p = .008) all agree, so the conclusion does not depend on the assumptions. That agreement is the strongest assurance you can report.

Remedies. Use Welch's t-test or Welch ANOVA by default (R oneway.test, Python anova_oneway(..., use_var="unequal")), with Games-Howell post hoc tests. If outliers or skew are severe and n is small, use a rank test (Mann-Whitney, Kruskal-Wallis) or a bootstrap CI for the mean difference (academic/stats-resampling-methods). For right-skewed positive outcomes (times, incomes) consider a log transformation (academic/stats-data-transformations).

Family 2: linear regression

The assumptions, remembered by the mnemonic LINE plus two: Linearity (the mean of y is a straight-line function of the predictors as specified), Independence of errors, Normality of errors (for small-sample inference and prediction intervals), Equal variance of errors (homoscedasticity), plus no undue influence from single observations and no severe collinearity between predictors.

The diagnostics.

fit_lm <- lm(score ~ hours + format, data = d)
summary(fit_lm)$coefficients
shapiro.test(residuals(fit_lm))
car::ncvTest(fit_lm)
car::vif(fit_lm)
diag_tab <- data.frame(id = d$id, rstudent = round(rstudent(fit_lm), 2),
                       leverage = round(hatvalues(fit_lm), 3), cooks = round(cooks.distance(fit_lm), 3))
head(diag_tab[order(-abs(diag_tab$rstudent)), ], 4)
car::outlierTest(fit_lm)
fit_lm2 <- update(fit_lm, subset = id != 24)
round(cbind(all = coef(fit_lm), without24 = coef(fit_lm2)), 2)
# par(mfrow = c(2, 2)); plot(fit_lm)       # the four standard diagnostic plots

(output abridged)
                Estimate Std. Error   t value     Pr(>|t|)
(Intercept)   40.7209172  5.3897812 7.5552079 2.788430e-07
hours          3.0961969  0.8832814 3.5053345 2.227439e-03
formatflipped  5.6705817  3.7975324 1.4932280 1.509863e-01
formatonline   0.8461969  3.4173689 0.2476165 8.069555e-01
Shapiro-Wilk normality test: W = 0.83329, p-value = 0.001087
Non-constant Variance Score Test (fitted values): Chisquare = 0.272451, Df = 1, p = 0.60169
           GVIF Df GVIF^(1/(2*Df))
hours  1.729306  1        1.315031
format 1.729306  2        1.146748
   id rstudent leverage cooks
24 24     5.11    0.129 0.430
10 10     1.30    0.159 0.077
20 20    -1.02    0.165 0.051
21 21    -1.00    0.237 0.077
   rstudent unadjusted p-value Bonferroni p
24 5.105406         6.2859e-05    0.0015086
                all without24
(Intercept)   40.72     41.91
hours          3.10      2.88
formatflipped  5.67      6.13
formatonline   0.85     -2.38

from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.outliers_influence import variance_inflation_factor, OLSInfluence
from statsmodels.stats.stattools import durbin_watson

fit_lm = smf.ols("score ~ hours + format", data=d).fit()
print(stats.shapiro(fit_lm.resid))
bp = het_breuschpagan(fit_lm.resid, fit_lm.model.exog)
print("Breusch-Pagan LM =", round(bp[0], 3), "p =", round(bp[1], 4))
X = fit_lm.model.exog
print("VIF:", [round(variance_inflation_factor(X, i), 3) for i in range(1, X.shape[1])])
print("Durbin-Watson =", round(durbin_watson(fit_lm.resid), 3))
inf = OLSInfluence(fit_lm)
tab = pd.DataFrame({"id": d["id"], "rstudent": inf.resid_studentized_external.round(2),
                    "leverage": inf.hat_matrix_diag.round(3), "cooks": inf.cooks_distance[0].round(3)})
print(tab.reindex(tab["rstudent"].abs().sort_values(ascending=False).index).head(4).to_string(index=False))

ShapiroResult(statistic=np.float64(0.8332853340039329), pvalue=np.float64(0.0010868978413425532))
Breusch-Pagan LM = 2.634 p = 0.4516
VIF: [np.float64(1.764), np.float64(1.429), np.float64(1.729)]
Durbin-Watson = 1.658
 id  rstudent  leverage  cooks
 24      5.11    0.129  0.430
 10      1.30    0.159  0.077
 20     -1.02    0.165  0.051
 21     -1.00    0.237  0.077

R and Python differ by design here, twice. First, R's car::ncvTest regresses the squared residuals on the fitted values (1 df: χ² = 0.27, p = .602), while statsmodels' het_breuschpagan with fit.model.exog uses all predictors (3 df: LM = 2.63, p = .452). The two functions also use different versions of the statistic: ncvTest is the original (non-studentised) Breusch-Pagan score test, and het_breuschpagan defaults to Koenker's studentised version (robust=True), which is less sensitive to non-normal residuals. Both say "no evidence of non-constant variance", but the numbers are not comparable. To reproduce R exactly in Python, use het_breuschpagan(fit_lm.resid, sm.add_constant(fit_lm.fittedvalues), robust=False) (χ² = 0.272, p = .602). Second, car::vif reports one generalised VIF (GVIF) for the whole format factor (1.73), while statsmodels gives one VIF per dummy column (1.764 for flipped, 1.429 for online) and 1.729 for hours. All are far below the usual concern thresholds of 5 to 10.

How to judge is shown step by step in Worked example A below, which walks through these residual diagnostics.

Remedies. For curvature, add a polynomial or spline term, or transform the predictor. For non-constant variance, use heteroscedasticity-consistent (HC3) standard errors (sandwich::vcovHC in R, fit = smf.ols(...).fit(cov_type="HC3") in Python) or transform the outcome. For influential points, never delete silently: check the record for errors, then report the model with and without the point. For collinearity, centre interaction terms, drop or combine redundant predictors, or accept wider CIs for those coefficients.

Family 3: generalised linear models

The assumptions. Independent observations; the link function and linear predictor are right (log odds linear in the predictors for logistic, log rate linear for Poisson); the variance follows the family (Poisson: variance equals the mean; binomial: variance p(1 − p)); and, for logistic regression, no complete separation and enough events per estimated parameter (about 10 is a common rule of thumb, more for stable estimates).

The overdispersion check (Poisson). Under a correct Poisson model, the Pearson chi-square statistic divided by the residual df should be close to 1:

$\hat{\phi} = \frac{1}{n - p} \sum_{i=1}^{n} \frac{(y_i - \hat{\mu}_i)^2}{\hat{\mu}_i}$

Values clearly above 1 (say above 1.5 with moderate n) mean the counts vary more than Poisson allows; the standard errors are then too small by a factor of about $\sqrt{\hat{\phi}}$.

fit_pois <- glm(visits ~ hours, family = poisson, data = d)
round(summary(fit_pois)$coefficients, 4)
disp <- sum(residuals(fit_pois, type = "pearson")^2) / df.residual(fit_pois)
cat("dispersion =", round(disp, 3), " df =", df.residual(fit_pois), "\n")
fit_nb <- MASS::glm.nb(visits ~ hours, data = d)
round(summary(fit_nb)$coefficients, 4); cat("theta =", round(fit_nb$theta, 3), "\n")
AIC(fit_pois, fit_nb)

fit_log <- glm(pass ~ hours, family = binomial, data = d)
round(summary(fit_log)$coefficients, 4)
table(pass = d$pass, hours7plus = d$hours >= 7)    # overlap check against separation

            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.2025     0.3458 -3.4776    5e-04
hours         0.4481     0.0441 10.1708    0e+00
dispersion = 2.79  df = 22 
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.8838     0.6280 -2.9998   0.0027
hours         0.5497     0.0913  6.0215   0.0000
theta = 2.445 
         df      AIC
fit_pois  2 134.8610
fit_nb    3 122.1455
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -6.4150     2.5947 -2.4724   0.0134
hours         1.0511     0.4261  2.4668   0.0136
    hours7plus
pass FALSE TRUE
   0    11    2
   1     4    7

fit_pois = smf.glm("visits ~ hours", data=d, family=sm.families.Poisson()).fit()
print(fit_pois.params.round(4))
print("dispersion =", round(fit_pois.pearson_chi2 / fit_pois.df_resid, 3))
fit_nb = smf.negativebinomial("visits ~ hours", data=d).fit(disp=0)
print(fit_nb.params.round(4), "theta =", round(1 / fit_nb.params["alpha"], 3))
fit_log = smf.logit("pass_ ~ hours", data=d).fit(disp=0)
print(fit_log.params.round(4), fit_log.pvalues.round(4))

Intercept   -1.2025
hours        0.4481
dtype: float64
dispersion = 2.79
Intercept   -1.8839
hours        0.5497
alpha        0.4090
dtype: float64 theta = 2.445
Intercept   -6.4150
hours        1.0511
dtype: float64 Intercept    0.0134
hours        0.0136
dtype: float64

How to judge. Library visits range from 0 to 25 with mean 6.25 and variance 47.5, so the variance is more than seven times the mean before any predictor is used. After fitting hours, the dispersion ratio is 2.79: the Poisson standard errors are too small by a factor of about $\sqrt{2.79} \approx 1.67$. The negative binomial model fits much better (AIC 122.1 against 134.9). Note the parameterisation difference: R's glm.nb reports θ = 2.445, statsmodels reports α = 0.409, and α = 1/θ. For the logistic model, passing and failing students overlap in study hours (some students with 7 or more hours failed, some with fewer passed), so there is no separation. With 11 passes and 13 fails there are about 5.5 events per parameter (two parameters), below the rule of thumb: the estimate is usable but imprecise, and adding more predictors would be unwise.

Remedies. Overdispersion: negative binomial (MASS::glm.nb, smf.negativebinomial) or quasi-Poisson (same estimates, SE inflated by $\sqrt{\hat{\phi}}$). Many zeros: consider a zero-inflated or hurdle model as a sensitivity analysis. Separation: penalised (Firth) logistic regression or exact methods. Few events: fewer predictors, or report the analysis as exploratory. Details in academic/stats-count-regression and academic/stats-logistic-regression.

Family 4: chi-square tests

The assumptions. Each unit appears in exactly one cell (independent units, no person counted twice), and the expected counts are large enough for the chi-square approximation. A common rule: all expected counts at least 5 for a 2 x 2 table; for larger tables, at least 80% of expected counts at least 5 and none below 1. Expected count for a cell is (row total × column total) / grand total.

tab <- table(d$format, d$pass)
tab
ct <- chisq.test(tab)
round(ct$expected, 2)
ct
fisher.test(tab)

(output abridged)
          0 1
  lecture 5 3
  flipped 1 7
  online  7 1
             0    1
  lecture 4.33 3.67
  flipped 4.33 3.67
  online  4.33 3.67
Pearson's Chi-squared test: X-squared = 9.3986, df = 2, p-value = 0.009102
Warning message: Chi-squared approximation may be incorrect
Fisher's Exact Test for Count Data: p-value = 0.0144

ct = pd.crosstab(d["format"], d["pass_"])
chi2, p, dof, exp = stats.chi2_contingency(ct)
print(exp.round(2)); print(round(chi2, 4), dof, round(p, 4))
mc = stats.MonteCarloMethod(n_resamples=99999, rng=np.random.default_rng(59))
print(stats.fisher_exact(ct.to_numpy(), method=mc).pvalue)

[[4.33 3.67]
 [4.33 3.67]
 [4.33 3.67]]
9.3986 2 0.0091
0.01474

How to judge. Every expected count is below 5 (4.33 and 3.67), so the chi-square p-value (.009) is not trustworthy; R even warns. Fisher's exact test gives p = .014. Here the conclusion survives, but the exact p-value is noticeably larger. R and Python differ by design: for tables larger than 2 x 2, R's fisher.test computes the exact p-value, while SciPy 1.18 estimates it by Monte Carlo simulation (0.0147 with this seed); without a seeded method the SciPy value changes slightly on every run.

Remedies. Fisher's exact test (or a Monte Carlo p-value); merge sparse categories only if the merge makes substantive sense and was planned; collect more data. If units appear more than once (same students before and after), use McNemar or Cochran's Q instead (academic/stats-contingency-tables).

Family 5: correlation

The assumptions. Pairs are independent; for Pearson's r, the relation is linear and no single point dominates (and for the usual CI, approximate bivariate normality); for Spearman's ρ, the relation is monotonic.

cor.test(d$hours, d$score)
cor.test(d$hours, d$score, method = "spearman", exact = FALSE)
cor.test(d$hours[-24], d$score[-24])      # sensitivity: without student 24

(output abridged)
Pearson:  t = 5.655, df = 22, p-value = 1.096e-05, r = 0.7696975, 95% CI 0.5312506 0.8951546
Spearman: S = 579.1, p-value = 2.62e-05, rho = 0.748216
Pearson without student 24: t = 7.758, df = 21, p-value = 1.343e-07, r = 0.8610098, 95% CI 0.6957289 0.9397022

r = stats.pearsonr(d["hours"], d["score"])
print(round(r.statistic, 4), r.pvalue, r.confidence_interval())
print(stats.spearmanr(d["hours"], d["score"]))

0.7697 1.0962953671991434e-05 ConfidenceInterval(low=np.float64(0.5312505924441787), high=np.float64(0.8951545758836629))
SignificanceResult(statistic=np.float64(0.7482160343994652), pvalue=np.float64(2.6201775238562153e-05))

How to judge. Pearson's r = .77, 95% CI [.53, .90]; Spearman's ρ = .75. Removing student 24 raises r to .86: one student out of 24 moves the correlation by .09. Report both, and say why student 24 was kept (a real, verified score) or excluded (a data error). Spearman's ρ is less affected by that point because it uses ranks.

Remedies. Spearman or Kendall for monotonic but curved relations or outliers; a bootstrap CI for r when the bivariate normality assumption is doubtful; and never infer causation from a correlation (academic/stats-correlation).

Family 6: repeated measures and mixed models

The assumptions. For repeated-measures ANOVA: complete data for every unit, normal residuals, and sphericity (the variances of all pairwise differences between conditions are equal). For linear mixed models: the random-effects structure matches the design (each student has their own intercept, and possibly their own slope), and both the residuals and the random effects are approximately normal. Mixed models do not assume sphericity and handle missing occasions.

mlm <- lm(cbind(q1, q2, q3) ~ 1, data = d)
mauchly.test(mlm, X = ~1)
anova(mlm, X = ~1, test = "Spherical")          # gives Greenhouse-Geisser and Huynh-Feldt
long <- reshape(d[, c("id", "q1", "q2", "q3")], direction = "long", varying = c("q1", "q2", "q3"),
                v.names = "quiz", timevar = "week", times = c(2, 6, 10), idvar = "id")
fit_mm <- lme4::lmer(quiz ~ factor(week) + (1 | id), data = long)
lme4::VarCorr(fit_mm)
shapiro.test(residuals(fit_mm))
shapiro.test(lme4::ranef(fit_mm)$id[, 1])

(output abridged)
Mauchly's test of sphericity: W = 0.73756, p-value = 0.03514
Greenhouse-Geisser epsilon: 0.7921
Huynh-Feldt epsilon:        0.8410
            Df      F num Df den Df     Pr(>F)     G-G Pr     H-F Pr
(Intercept)  1 46.407      2     46 9.2792e-12 9.1555e-10 3.1048e-10
 Groups   Name        Std.Dev.
 id       (Intercept) 1.08152 
 Residual             0.55441 
Shapiro-Wilk (residuals): W = 0.98026, p-value = 0.3178
Shapiro-Wilk (random intercepts): W = 0.97527, p-value = 0.7955

Y = d[["q1", "q2", "q3"]].to_numpy(float)
K = np.array([[1, -1, 0], [1, 1, -2]], float)
K = (K.T / np.linalg.norm(K, axis=1)).T          # orthonormal contrasts
S = np.cov(Y @ K.T, rowvar=False)
k, n = 2, Y.shape[0]
W = np.linalg.det(S) / (np.trace(S) / k) ** k
f = 1 - (2 * k**2 + k + 2) / (6 * k * (n - 1))
chi = -(n - 1) * f * np.log(W); df = k * (k + 1) / 2 - 1
gg = np.trace(S) ** 2 / (k * np.trace(S @ S))
print("Mauchly W =", round(W, 4), "p =", round(stats.chi2.sf(chi, df), 4), "GG epsilon =", round(gg, 4))
long = d.melt(id_vars="id", value_vars=["q1", "q2", "q3"], var_name="week", value_name="quiz")
fit_mm = smf.mixedlm("quiz ~ C(week)", long, groups=long["id"]).fit(reml=True)
print("random intercept SD =", round(float(np.sqrt(fit_mm.cov_re.iloc[0, 0])), 4),
      "residual SD =", round(float(np.sqrt(fit_mm.scale)), 4))
print(stats.shapiro(fit_mm.resid))

Mauchly W = 0.7376 p = 0.0351 GG epsilon = 0.7921
random intercept SD = 1.0815 residual SD = 0.5544
ShapiroResult(statistic=np.float64(0.9802639273299788), pvalue=np.float64(0.31786711122268063))

In the R table, the row labelled (Intercept) is the test of the quiz-occasion effect (the within-student contrasts are tested against zero), with uncorrected, Greenhouse-Geisser and Huynh-Feldt p-values. statsmodels has no built-in Mauchly test (its AnovaRM reports uncorrected F only), so the Python block computes Mauchly's W and the Greenhouse-Geisser ε from orthonormal contrasts; the results match R exactly.

How to judge. Mauchly's test is significant (W = 0.74, p = .035) and Greenhouse-Geisser ε = 0.79, so sphericity is doubtful. Like normality tests, Mauchly's test has low power in small samples and too much power in large ones, so many analysts apply the Greenhouse-Geisser correction routinely. Here the quiz effect is huge either way (corrected p < .001). In the mixed model, students differ a lot in their overall level (random intercept SD 1.08 points) relative to the week-to-week noise (residual SD 0.55), and neither the residuals nor the random intercepts show a normality problem.

Remedies. Report Greenhouse-Geisser-corrected df and p-values (or Huynh-Feldt when ε is above about 0.75); move to a linear mixed model when data are missing or timing varies; add random slopes when students differ in their rate of change (academic/stats-repeated-measures, academic/stats-mixed-models).

Worked example A (fully solved): a residual diagnostics walk-through

Question. In the model score ~ hours + format, are the assumptions good enough to report the hours coefficient (3.10 points per weekly study hour, p = .002)?

Step 1: residuals vs fitted (linearity and equal variance). In a good plot the points form a horizontal band around zero with no curve and no funnel. The chart below shows that 23 students sit in a band from about −6 to +8 points, with no curve and no funnel. Student 24 sits alone at +20.95. Decision: linearity is fine; the band is even; one outlier needs attention.

Step 2: QQ plot (normality). In a good QQ plot the points follow the diagonal line. Here 23 points do; student 24 is far above the line at the top right. Shapiro-Wilk agrees (W = 0.833, p = .001), but the plot tells you why: one point, not a skewed distribution. Decision: the non-normality is caused by a single outlier, so treat it as an outlier problem, not a transformation problem. Without student 24, Shapiro-Wilk gives W = 0.956, p = .390.

Step 3: scale-location (equal variance). The square root of the absolute standardised residuals against fitted values should show a flat trend. The score test agrees there is no funnel (χ²(1) = 0.27, p = .602). Decision: equal variance is acceptable.

Step 4: influence (residuals vs leverage, Cook's distance). Three numbers per point: the studentised residual (how surprising the point is, roughly on a t scale), the leverage (how unusual its predictor values are; the average is p/n = 4/24 = 0.167) and Cook's distance (how much all fitted values change if the point is removed). Student 24 has a studentised residual of 5.11 (Bonferroni-adjusted p = .002 in car::outlierTest), ordinary leverage 0.129 and Cook's distance 0.430. Cook's distance is above the screening cut-off of 4/n = 0.167 but below 1. The formula is

$D_i = \frac{e_i^2}{p \, s^2} \cdot \frac{h_{ii}}{(1 - h_{ii})^2}$

where $e_i$ is the residual, $p$ the number of coefficients (4), $s^2$ the residual variance (6.602² = 43.59) and $h_{ii}$ the leverage. By hand: $20.95^2 / (4 \times 43.59) = 438.9 / 174.4 = 2.517$ and $0.129 / 0.871^2 = 0.170$, so $D_{24} \approx 2.517 \times 0.170 = 0.43$, matching the output. Student 24 is an outlier in y with ordinary x values: influential on the online coefficient, not on the slope.

Step 5: sensitivity analysis and decision. Refit without student 24 and compare: the hours coefficient moves from 3.10 to 2.88 (a change of 0.22 points per hour); the online vs lecture difference moves from +0.85 to −2.38, a change of sign, though neither is statistically distinguishable from zero. Decision: (1) check the record for student 24; if the score is a data-entry error, correct it; (2) if the score is real, keep it, report the main model, and report the model without student 24 as a sensitivity analysis; (3) the conclusion about study hours is robust, the conclusion about format differences is fragile.

Report sentence. "Residual diagnostics showed one outlying student (studentised residual 5.11, Cook's D = 0.43). Excluding this student changed the study-hours coefficient only slightly (3.10 to 2.88 points per hour) but reversed the sign of the small online vs lecture difference; both models are reported."

Notation used above:

Worked example B (half solved): is the Poisson model good enough?

Question. Can the Poisson model visits ~ hours be used to report the rate ratio for study hours?

1. Before fitting, compare the mean and variance of visits: mean 6.25, variance 47.5. For a Poisson variable these should be roughly [ ].
2. After fitting, the dispersion ratio is 61.38 / 22 = [ ]. This is well above 1, so the Poisson standard errors are [ ] (too small / too large) by a factor of about [ ].
3. The negative binomial model has AIC 122.1 against 134.9 for Poisson, so [ ] fits better.
4. The rate ratio per extra study hour from the negative binomial model is exp(0.5497) = 1.733, 95% Wald CI [1.449, 2.072]. In words: [ ].
5. Decision: report the [ ] model, and mention the Poisson result as [ ].

Exercise C (unsolved): expected counts by hand

Harbour City Library (fictional) records the loan type of 60 loans at two branches.

(a) Compute all six expected counts. (b) Is the chi-square test of independence acceptable? (c) Compute the chi-square statistic and its df. (d) Compute Cramér's V. (e) What would you do if one branch had only 10 loans?

Common mistakes and misconceptions

Self-test questions

1. In Family 1, Shapiro-Wilk on the ANOVA residuals gave p = .234 but the regression residuals gave p = .001 for the same students. Why can both be true?
2. Rank these from most to least serious when violated in a regression with n = 200: non-normal residuals, students clustered in 10 classes analysed as independent, mild non-constant variance.
3. A Poisson regression has a dispersion ratio of 3.2. What does that do to the reported p-values, and name two remedies.
4. Why does SciPy's fisher_exact give a slightly different p-value each run for a 3 x 2 table, and how do you make it reproducible?
5. Which two numbers would you report to show that one observation is influential, and what is the screening cut-off for Cook's distance used here?
6. Mauchly's test gives p = .035 and Greenhouse-Geisser ε = 0.79. What do you report?

My checks

My notes and next steps

Sources and further reading

Based on the NIST/SEMATECH e-Handbook of Statistical Methods (sections on assumptions, residual analysis and goodness of fit) and the original papers below, written in original words and examples.

• NIST/SEMATECH e-Handbook of Statistical Methods, section 1.3.5.10 Levene Test and section 4.4.4 How can I tell if a model fits my data?, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain.
• Box, G. E. P., and Cox, D. R. (1964). An analysis of transformations. Journal of the Royal Statistical Society, Series B, 26(2), 211–243.
• Cook, R. D. (1977). Detection of influential observation in linear regression. Technometrics, 19(1), 15–18.
• Levene, H. (1960). Robust tests for equality of variances. In I. Olkin et al. (Eds.), Contributions to Probability and Statistics: Essays in Honor of Harold Hotelling (pp. 278–292). Stanford University Press.
• Brown, M. B., and Forsythe, A. B. (1974). Robust tests for the equality of variances. Journal of the American Statistical Association, 69(346), 364–367.
• Shapiro, S. S., and Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3–4), 591–611.
• Breusch, T. S., and Pagan, A. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5), 1287–1294.
• Mauchly, J. W. (1940). Significance test for sphericity of a normal n-variate distribution. Annals of Mathematical Statistics, 11(2), 204–209.
• Greenhouse, S. W., and Geisser, S. (1959). On methods in the analysis of profile data. Psychometrika, 24(2), 95–112.
• Nelder, J. A., and Wedderburn, R. W. M. (1972). Generalized linear models. Journal of the Royal Statistical Society, Series A, 135(3), 370–384.
• Welch, B. L. (1947). The generalization of "Student's" problem when several different population variances are involved. Biometrika, 34(1–2), 28–35.
• Fox, J., and Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). Sage. The companion book of the car package (ncvTest, vif with GVIF, outlierTest, leveneTest).
• Further reading (by name only): de Smith, M. J., Statistical Analysis Handbook (2024 edition); Diez, Çetinkaya-Rundel and Barr, OpenIntro Statistics.

Answer key

Example B.

1. Equal (for a Poisson variable the variance equals the mean).
2. 61.38 / 22 = 2.79. The Poisson standard errors are too small by a factor of about $\sqrt{2.79} \approx$ 1.67.
3. The negative binomial model fits better (lower AIC by about 12.7).
4. Each extra weekly study hour is associated with about 73% more library visits (rate ratio 1.73, 95% CI [1.45, 2.07]).
5. Report the negative binomial model; mention the Poisson result as a comparison or sensitivity analysis, noting its overdispersion (dispersion ratio 2.79).

Exercise C.

• (a) North and South have the same totals (30), so both rows have the same expected counts: fiction 30 × 30 / 60 = 15, non-fiction 30 × 20 / 60 = 10, children's 30 × 10 / 60 = 5.
• (b) Yes: the smallest expected count is exactly 5, so every cell meets the rule.
• (c) χ² = (18 − 15)²/15 + (9 − 10)²/10 + (3 − 5)²/5 + (12 − 15)²/15 + (11 − 10)²/10 + (7 − 5)²/5 = 0.6 + 0.1 + 0.8 + 0.6 + 0.1 + 0.8 = 3.0, df = (2 − 1)(3 − 1) = 2, p = .223 (chisq.test gives X-squared = 3, df = 2, p-value = 0.2231).
• (d) Cramér's V = $\sqrt{3.0 / (60 \times 1)}$ = 0.22 (the 1 is min(rows − 1, columns − 1)).
• (e) The children's expected counts would drop below 5 for that branch: use Fisher's exact test (or a Monte Carlo p-value), or merge categories only if that answers the question.

Self-test.

1. The ANOVA residuals measure each score from its group mean, and student 24's high score inflates the online group mean and SD, partly hiding the outlier. The regression also uses hours, which predicts most students well, so student 24's residual (20.95) stands out sharply against a residual SD of 6.6. Different models, different residuals.
2. Most serious: clustering ignored (SE too small, no sample size fixes it); then mild non-constant variance (use HC3 SE); least: non-normal residuals with n = 200 (the central limit theorem protects the coefficient tests).
3. Standard errors are too small by about $\sqrt{3.2} \approx 1.79$, so p-values are too small and CIs too narrow. Remedies: negative binomial regression; quasi-Poisson; also check for a missing predictor or excess zeros.
4. For tables larger than 2 x 2, SciPy estimates the p-value by Monte Carlo simulation. Pass a seeded method, e.g. method=stats.MonteCarloMethod(n_resamples=99999, rng=np.random.default_rng(59)), or use R's exact fisher.test.
5. The studentised residual and Cook's distance (with leverage as context); the screening cut-off used here is 4/n (0.167 for n = 24), with values near or above 1 as clearly influential.
6. Report the Greenhouse-Geisser corrected test: multiply both df by ε (2 × 0.79 = 1.58 and 46 × 0.79 = 36.4), and give the corrected p-value with ε, e.g. "F(1.58, 36.4) = 46.41, p < .001, ε = 0.79".

Revision history