# 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 | Field | Value | |---|---| | Course or project | [STAT 201 / thesis chapter 4 / client project] | | Learner or analyst | [Your name] | | Data set | [Lakemoor University study-habits data, 24 students (fictional) / your data set and source] | | Software and versions | R 4.5.1 (stats) with car 3.1.5, lme4 2.1.0, MASS 7.3.65; Python 3.14 with numpy 2.5.4, pandas 3.1.0, scipy 1.18.1, statsmodels 0.15.0 | | Linked analysis plan | [academic/stats-analysis-plan, version 0.x] | | Date | YYYY-MM-DD | | Version | [0.1] | ## 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. ```r 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") ``` ```text n = 15, skewed (lognormal): rejected 0.508 n = 5000, near-normal (t, df = 10): rejected 1 ``` ```python 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) ``` ```text 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. | id | format | hours | score | visits | q1 | q2 | q3 | |---|---|---|---|---|---|---|---| | 1–8 | lecture | 4, 6, 5, 8, 3, 7, 5, 6 | 56, 55, 57, 66, 48, 57, 62, 61 | 0, 9, 1, 15, 0, 2, 6, 4 | 5, 6, 5, 7, 4, 6, 5, 6 | 6, 6, 7, 8, 5, 6, 6, 7 | 7, 8, 6, 8, 6, 8, 7, 7 | | 9–16 | flipped | 7, 9, 6, 10, 8, 5, 9, 7 | 62, 82, 68, 76, 75, 57, 69, 71 | 3, 18, 1, 25, 12, 0, 8, 10 | 6, 7, 5, 8, 6, 5, 7, 6 | 6, 8, 7, 9, 8, 6, 7, 7 | 8, 9, 7, 9, 8, 7, 9, 8 | | 17–24 | online | 3, 5, 4, 6, 2, 7, 4, 5 | 50, 51, 58, 54, 42, 58, 53, 78 | 0, 7, 0, 11, 0, 15, 1, 2 | 4, 5, 5, 6, 3, 6, 4, 5 | 5, 4, 6, 7, 4, 5, 5, 6 | 5, 6, 6, 6, 5, 7, 5, 7 | ```r # 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) ``` ```python 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 | Assumptions (most important first) | Main plots | Supporting tests | First remedy | |---|---|---|---|---| | t-tests and ANOVA | Independence; no extreme outliers; approx. normal residuals; equal variances (classic F and pooled t only) | Dot plot or box plot per group; QQ plot of residuals | Levene (median-centred); Shapiro-Wilk on residuals | Welch t or Welch ANOVA; rank test; bootstrap | | Linear regression | Independence; linearity; constant variance; no undue influence; normal residuals; no severe collinearity | Residuals vs fitted; QQ; scale-location; residuals vs leverage with Cook's distance | Breusch-Pagan; Durbin-Watson; VIF | Transform or add terms; robust (HC3) SE; report with and without influential points | | GLMs (logistic, Poisson) | Independence; correct link and linear predictor; variance function (no overdispersion); no separation; enough events | Binned residuals; observed vs fitted counts | Pearson dispersion ratio; events per parameter | Negative binomial or quasi-Poisson; penalised logistic; fewer predictors | | Chi-square tests | Independent units (each counted once); expected counts large enough | The expected-count table itself | None needed: compute expected counts | Fisher's exact; merge sparse categories with a reason; McNemar for paired data | | Correlation | Linear (Pearson) or monotonic (Spearman) relation; no influential outliers; independent pairs | Scatter plot | None needed | Spearman or Kendall; report with and without the outlier | | Repeated measures and mixed models | Correct grouping structure; sphericity (RM ANOVA); normal residuals and random effects | Residuals vs fitted; QQ of residuals and of random effects; spaghetti plot | Mauchly's test; Greenhouse-Geisser ε | Greenhouse-Geisser correction; mixed model; random slopes | ```mermaid flowchart TD A["Choose method (academic/stats-test-selector)"] --> B["List its assumptions from the overview table"] B --> C{"Independence and design correct?"} C -->|No| C1["Change the model: paired test, mixed model, cluster-robust SE"] C -->|Yes| D["Fit the model, then plot residuals"] D --> E{"Departure large enough to change the conclusion?"} E -->|No| F["Proceed and report the checks"] E -->|Yes| G["Apply the planned remedy"] G --> H["Report main and sensitivity results"] F --> H ``` ## 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. ```r 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 ``` ```text (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 ``` ```python 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)) ``` ```text (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: **L**inearity (the mean of y is a straight-line function of the predictors as specified), **I**ndependence of errors, **N**ormality of errors (for small-sample inference and prediction intervals), **E**qual variance of errors (homoscedasticity), plus no undue influence from single observations and no severe collinearity between predictors. **The diagnostics.** | Assumption | Plot | Supporting number | R | Python | |---|---|---|---|---| | Linearity | Residuals vs fitted; residuals vs each predictor | None | `plot(fit, which = 1)` | `plt.scatter(fit.fittedvalues, fit.resid)` | | Normality | QQ plot of residuals | Shapiro-Wilk | `plot(fit, which = 2)`, `shapiro.test(residuals(fit))` | `sm.qqplot(fit.resid, line="s")`, `stats.shapiro(fit.resid)` | | Equal variance | Scale-location plot | Breusch-Pagan / score test | `plot(fit, which = 3)`, `car::ncvTest(fit)` | `het_breuschpagan(fit.resid, fit.model.exog)` | | Independence | Residuals in time or collection order | Durbin-Watson (only for ordered data) | `car::durbinWatsonTest(fit)` | `durbin_watson(fit.resid)` | | Influence | Residuals vs leverage with Cook's distance | Cook's D, studentised residuals, leverage | `plot(fit, which = 5)`, `cooks.distance(fit)` | `OLSInfluence(fit).cooks_distance` | | Collinearity | None | VIF (GVIF for factors) | `car::vif(fit)` | `variance_inflation_factor(X, i)` | ```r 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 ``` ```text (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 ``` ```python 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)) ``` ```text 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}}$. ```r 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 ``` ```text 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 ``` ```python 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)) ``` ```text 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. ```r tab <- table(d$format, d$pass) tab ct <- chisq.test(tab) round(ct$expected, 2) ct fisher.test(tab) ``` ```text (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 ``` ```python 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) ``` ```text [[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. ```r 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 ``` ```text (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 ``` ```python r = stats.pearsonr(d["hours"], d["score"]) print(round(r.statistic, 4), r.pvalue, r.confidence_interval()) print(stats.spearmanr(d["hours"], d["score"])) ``` ```text 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. ```r 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]) ``` ```text (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 ``` ```python 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)) ``` ```text 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. ```vega-lite { "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "description": "Residuals against fitted values for the model score ~ hours + format; student 24 stands out", "width": 360, "height": 240, "data": {"values": [ {"id": 1, "fitted": 53.11, "resid": 2.89}, {"id": 2, "fitted": 59.3, "resid": -4.3}, {"id": 3, "fitted": 56.2, "resid": 0.8}, {"id": 4, "fitted": 65.49, "resid": 0.51}, {"id": 5, "fitted": 50.01, "resid": -2.01}, {"id": 6, "fitted": 62.39, "resid": -5.39}, {"id": 7, "fitted": 56.2, "resid": 5.8}, {"id": 8, "fitted": 59.3, "resid": 1.7}, {"id": 9, "fitted": 68.06, "resid": -6.06}, {"id": 10, "fitted": 74.26, "resid": 7.74}, {"id": 11, "fitted": 64.97, "resid": 3.03}, {"id": 12, "fitted": 77.35, "resid": -1.35}, {"id": 13, "fitted": 71.16, "resid": 3.84}, {"id": 14, "fitted": 61.87, "resid": -4.87}, {"id": 15, "fitted": 74.26, "resid": -5.26}, {"id": 16, "fitted": 68.06, "resid": 2.94}, {"id": 17, "fitted": 50.86, "resid": -0.86}, {"id": 18, "fitted": 57.05, "resid": -6.05}, {"id": 19, "fitted": 53.95, "resid": 4.05}, {"id": 20, "fitted": 60.14, "resid": -6.14}, {"id": 21, "fitted": 47.76, "resid": -5.76}, {"id": 22, "fitted": 63.24, "resid": -5.24}, {"id": 23, "fitted": 53.95, "resid": -0.95}, {"id": 24, "fitted": 57.05, "resid": 20.95} ]}, "layer": [ { "mark": {"type": "rule", "strokeDash": [4, 4], "color": "#888888"}, "encoding": {"y": {"datum": 0}} }, { "mark": {"type": "point", "filled": true, "size": 60}, "encoding": { "x": {"field": "fitted", "type": "quantitative", "title": "Fitted exam score", "scale": {"zero": false}}, "y": {"field": "resid", "type": "quantitative", "title": "Residual (points)"}, "color": {"condition": {"test": "datum.id === 24", "value": "#d62728"}, "value": "#4c78a8"}, "tooltip": [{"field": "id"}, {"field": "fitted"}, {"field": "resid"}] } } ] } ``` **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: | Symbol | Meaning | |---|---| | $e_i$ | Residual for observation i: observed minus fitted | | $h_{ii}$ | Leverage of observation i (diagonal of the hat matrix) | | $s$ | Residual standard error (6.602 here) | | $p$ | Number of estimated coefficients including the intercept | | $D_i$ | Cook's distance of observation i | | $\hat{\phi}$ | Dispersion ratio for a GLM | ## 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. | Branch | Fiction | Non-fiction | Children's | Total | |---|---|---|---|---| | North | 18 | 9 | 3 | 30 | | South | 12 | 11 | 7 | 30 | | Total | 30 | 20 | 10 | 60 | (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 | Wrong belief | Correct statement | |---|---| | "The outcome must be normally distributed." | For t-tests, ANOVA and regression, the residuals (and, more importantly, the sampling distribution of the estimates) should be approximately normal. The raw outcome pooled across groups can look non-normal when groups differ. | | "Shapiro-Wilk p > .05 proves normality." | A non-significant test means the test could not detect a departure, which is common with small n. It is not proof. | | "Shapiro-Wilk p < .05 means I must use a nonparametric test." | With large n, trivial departures are significant; the t-test and regression are robust to them. Look at the plot and the size of the departure. | | "Levene's test is not significant, so variances are equal." | With small groups Levene has low power. Use Welch by default and the test does not matter. | | "Outliers should be removed." | Outliers should be investigated. Remove only data errors; otherwise report results with and without them. | | "A GLM has no assumptions because it does not assume normality." | GLMs assume a mean-variance relation, a link function, independence and (for logistic) no separation. | | "If the chi-square observed counts are all above 5, the test is fine." | The rule is about **expected** counts. | | "Nonparametric tests are assumption-free." | They still assume independence, and Mann-Whitney compares medians only when the two distributions have the same shape. | | "Transforming the data always fixes non-normality." | A log transformation changes the question (ratios of geometric means, not differences of means) and cannot fix an outlier caused by a data error. | | "Durbin-Watson checks independence for any data." | It checks first-order autocorrelation in data with a meaningful order (time, sequence). For clustered data, the design tells you, not a test. | ## 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 | Method | Assumption | Check used (plot and number) | Result | Decision | Planned in advance? | |---|---|---|---|---|---| | *Example:* Linear regression | Influence | Residuals vs leverage; Cook's D | Student 24: D = 0.43, rstudent 5.11 | Keep, report with and without | Yes | | *Example:* Poisson regression | No overdispersion | Pearson dispersion ratio | 2.79 | Switch to negative binomial | Yes | | [Method] | [Assumption] | [Plot / statistic] | [Value] | [Decision] | [Yes / No, explain] | ## My notes and next steps | Item | My entry | |---|---| | Confidence reading diagnostic plots (1 to 5) | [ ] | | Family I find hardest to check | [ ] | | Questions for my supervisor or tutor | [ ] | | Mistakes logged in `academic/stats-error-log` | [link or count] | | Next review date | YYYY-MM-DD | ## 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 | Version | Date | Author | Change summary | |---|---|---|---| | 0.1 | YYYY-MM-DD | [Your name] | First draft from the template | | 0.2 | YYYY-MM-DD | [Your name] | *Example:* Replaced example data with thesis data; negative binomial chosen after dispersion check |