ANOVA in Python with statsmodels: anova_lm, Type III Contrasts, and the AnovaRM Limits

statsmodels has no ANOVA dialog box to keep you honest. It will accept almost any formula you hand it, fit something, and print a tidy F table — and there are three specific ways that table can be answering a different question than the one you asked. None of them raise an error. All three are visible in the output if you know which column to look at.

This page walks the working path for a one-way and a factorial ANOVA in statsmodels, then the repeated-measures case, where the library’s limits are hard enough that you may need a different package.

The three things that silently go wrong

1. A grouping column stored as integers, without C(). patsy treats it as a continuous predictor and fits a straight line through your group codes. You get an F test on 1 numerator degree of freedom instead of k − 1.
2. anova_lm defaults to Type I. The documented default for typ is I — sequential sums of squares. On an unbalanced factorial design the answer then depends on the order you wrote the terms in the formula.
3. typ=3 with patsy’s default contrasts. Type III main effects are only interpretable under sum-to-zero-style coding, and patsy’s default is Treatment (dummy) coding. This combination runs cleanly and reports the wrong thing.

Everything below is about avoiding those three.

The minimum correct one-way ANOVA

Use the formula API, wrap the factor in C(), and pass the fitted result to anova_lm:

• import statsmodels.api as sm
• from statsmodels.formula.api import ols
• model = ols('score ~ C(treatment)', data=df).fit()
• sm.stats.anova_lm(model, typ=2)

The C() wrapper is not decoration. patsy’s documentation describes C() as marking data as categorical “including data which would not automatically be treated as categorical, such as a column of integers”. If treatment holds 0/1/2 or dose values of 0/10/20, then score ~ treatment fits a regression slope, not an ANOVA.

The one-line check: look at the df column of the ANOVA table for your factor. It must equal the number of levels minus one. A factor row showing 1.0 when you have four groups means patsy read it as numeric.

Reading the table you get back

The columns differ by type, which trips people who write code against them. typ=1 returns df, sum_sq, mean_sq, F and PR(>F). typ=2 and typ=3 return sum_sq, df, F and PR(>F) — no mean_sq column at all. Type III output also carries an extra Intercept row, which is a test that the grand mean differs from zero and is almost never of interest.

For a genuine one-way design the SS type is irrelevant: with a single factor there is nothing for the sequential ordering to depend on, and Types I, II and III agree. pingouin’s documentation states the same thing for its ss_type argument, which “has no impact on one-way designs or balanced N-way ANOVA”. If your design is one-way and balanced, stop worrying about types and go read the F-statistic itself.

Factorial designs: which typ, and what typ=3 costs you

The moment you have two or more factors and unequal cell counts, the SS type changes the answer. The statsmodels documentation is explicit about each:

• Type I (sequential). “Each term is tested after the terms that precede it in the model. Consequently, the results depend on the order of the terms when the design is unbalanced.” This is the default, and it is a reporting hazard: swapping two terms in your formula string changes your p-values.
• Type II. “Each term is tested after all other terms except higher-order terms that contain it.” This respects the marginality principle.
• Type III. “Each term is tested after all other terms in the model, including higher-order terms that contain it” — and, critically, “may depend on the contrast coding used for categorical factors”.

The contrast requirement, spelled out

That last clause is the whole problem, and the R car package documents it more bluntly than statsmodels does. Type III tests “will normally only be sensible when using contrasts that, for different terms, are orthogonal in the row-basis of the model, such as those produced by contr.sum, contr.poly, or contr.helmert, but not by the default contr.treatment”. patsy’s default is Treatment coding — “simply giving Patsy a categorical variable causes it to be coded using the default Treatment coding scheme”.

So the two defaults compose into exactly the invalid combination: patsy gives you dummy coding, and typ=3 then reports what is effectively a simple effect at the reference level of the other factor, labelled as a main effect. The fix is in the formula, not the anova_lm call. statsmodels’ own ANOVA user guide writes it this way:

• moore_lm = ols('conformity ~ C(fcategory, Sum)*C(partner_status, Sum)', data=data).fit()
• table = sm.stats.anova_lm(moore_lm, typ=2)

Note what the official example does: it sets Sum coding and asks for Type II. Type II is the safer default of the two, because — per the same car documentation — “type-II tests are invariant with respect to (full-rank) contrast coding”. You cannot get a Type II test wrong by forgetting to change your contrasts.

Practical rule: use typ=2 unless you specifically need to test a main effect in the presence of an interaction that you intend to keep in the model. If you do need Type III, switch every categorical term to C(x, Sum) first, and say so when you report. The same logic, and the same trap, applies when running ANOVA in R — the packages differ, the statistics do not.

The robust argument that does nothing on the default

anova_lm accepts robust in {None, 'hc0', 'hc1', 'hc2', 'hc3'}, with hc3 recommended. But in the statsmodels source, robust is “accepted for interface consistency but not used for Type I sums of squares”. Passing robust='hc3' without also passing typ=2 or typ=3 is a no-op that returns no warning.

Unequal variances: statsmodels has two different one-way defaults

This catches people who mix the two APIs. anova_lm on an OLS fit assumes homoscedasticity, as does scipy.stats.f_oneway, whose documented assumptions are that “the samples are independent”, “each sample is from a normally distributed population” and “the population standard deviations of the groups are all equal”.

But statsmodels.stats.oneway.anova_oneway has the opposite default. Its signature is anova_oneway(data, groups=None, use_var='unequal', welch_correction=True, trim_frac=0), and use_var='unequal' — the default — “implements Welch Anova with Satterthwaite-Welch degrees of freedom”. To get the classical test you must pass use_var='equal', described as “the standard Anova”. A third option, use_var='bf', gives the Brown-Forsythe method with Mehrotra corrected degrees of freedom, and trim_frac switches to a trimmed-mean and winsorized-variance (Yuen) variant.

So “I ran a one-way ANOVA in statsmodels” is ambiguous between two functions that disagree by default. State which one you used. If variances are clearly unequal, the Welch default is the right one — and note that scipy.stats.f_oneway now also takes an equal_var argument, while the SciPy docs point to kruskal and alexandergovern as alternatives “although with some loss of power” when assumptions fail.

Post hoc comparisons

pairwise_tukeyhsd(endog, groups, alpha=0.05, use_var='equal') is the built-in. Two things about it are easy to miss.

First, it does not take your fitted model. It takes a response array and a single grouping array, so for a factorial design you have to build a combined cell-label column yourself — and once you do, you are comparing cells, not marginal means, which is a different question from the one the factorial F tests answered.

Second, its default is fragile in exactly the situation people reach for it. The documentation warns that “Tukey-hsd is not robust to heteroscedasticity, i.e., variance differ across groups, especially if group sizes also vary. In those cases, the actual size (rejection rate under the Null hypothesis) might be far from the nominal size of the test.” Passing use_var='unequal' switches it to the Games-Howell procedure, which uses Welch’s t-test for unequal variances and “approximately maintains size unless samples are very small”. That one keyword is the difference between a defensible post hoc and an inflated false-positive rate.

Repeated measures: know AnovaRM’s three hard limits before you start

AnovaRM(data, depvar, subject, within=None, between=None, aggregate_func=None) handles within-subject designs from long-format data. It has three constraints, and the third is the one that damages a manuscript.

1. Fully balanced designs only. The documentation states the implementation “currently only supports fully balanced designs”; the source raises ValueError with the message Data is unbalanced. One dropped subject-by-condition observation and it stops.
2. No between-subject factors. The between argument exists but is not implemented; supplying it raises NotImplementedError("Between subject effect not yet supported!"). A mixed within/between design cannot be run here at all.
3. No sphericity correction. The class docstring says plainly that “calculation of between-subject effects and corrections for violation of sphericity are not yet implemented”. The output table’s columns are F Value, Num DF, Den DF and Pr > F — there is no epsilon, no Greenhouse-Geisser corrected p-value, and no Mauchly test.

That third point is why a clean-looking AnovaRM table is not automatically reportable. If your within-subject factor has three or more levels, sphericity has to be checked and, if violated, corrected, and nothing in statsmodels will do it or warn you that it was skipped. Reviewers in psychology and neuroscience routinely ask for the epsilon.

One more operational detail: if your data hold more than one observation per subject per cell, they “need to be aggregated into a single observation”, and an exception is raised if aggregation is required but no aggregate_func was supplied. Passing aggregate_func='mean' silences that — make sure averaging is actually what you want, because it quietly discards within-cell variance.

When to switch to pingouin

pingouin.rm_anova(data, dv, within, subject, correction='auto', detailed=False, effsize='ng2') fills the gap. With correction='auto' it computes Mauchly’s test to decide whether p-values need correcting, and the returned frame includes W_spher, p_spher, sphericity, eps and p_GG_corr alongside the usual Source, ddof1, ddof2, F and p_unc. For the mixed design AnovaRM refuses, pingouin.mixed_anova(data, dv, within, subject, between, correction='auto', effsize='np2') takes one within and one between factor and reports the same sphericity columns.

The tradeoff to know about: pingouin uses “a strict listwise approach (= complete-case analysis)”, so any subject with a missing cell is removed entirely before the test, which the documentation warns “could drastically decrease the power of the ANOVA if many missing values are present”. Count your surviving subjects, do not assume your n is what the dataframe says.

Effect size: statsmodels will not give you one

There is no eta-squared column anywhere in anova_lm’s output. You compute it from the table you already have:

• Eta-squared is the effect sum of squares over the total sum of squares — in code, table['sum_sq'] / table['sum_sq'].sum(). Because all effects sum to the total, these are comparable within one study.
• Partial eta-squared is the effect sum of squares over the sum of the effect and error sums of squares — ss_effect / (ss_effect + ss_resid), using the Residual row.

Lakens (2013) is worth reading before you pick one. Partial eta-squared “differs when the same two means are compared in a within-subjects design or a between-subjects design”, so it is not comparable across studies with different designs; plain eta-squared is also design-dependent because “the total variability in a study depends on the design of a study, and increases when additional variables are manipulated”. Omega-squared has “been suggested to correct for” the upward bias in these estimates, though it is “at best a less biased estimate” and the difference shrinks with larger samples.

If you would rather not hand-roll it, pingouin.anova(data, dv, between, ss_type=2, detailed=False, effsize='np2') returns np2 (partial eta-squared) as a column by default, with effsize='n2' for plain eta-squared. Note pingouin’s ss_type defaults to 2, not 1 — another reason two “Python ANOVA” scripts can disagree.

What to put in the paper

Because the defaults differ across functions and across packages, a reader cannot reconstruct your F from the numbers alone. Report:

• The F ratio with both degrees of freedom, the p-value, and an effect size with its type named.
• Which sums-of-squares type you used. “Type II sums of squares” is one clause and removes the largest source of ambiguity.
• The contrast coding, if and only if you used Type III. If you used typ=3 without C(x, Sum), the honest fix is to rerun it, not to omit the detail.
• Whether the equal-variance assumption held, and which test you used if it did not (Welch anova_oneway, Brown-Forsythe, Kruskal-Wallis).
• For repeated measures, the sphericity test result and the correction applied — which means you did not get the numbers from AnovaRM alone.

Frequently asked questions

Why does my Python ANOVA disagree with SPSS?

The most common cause on an unbalanced factorial design is the sums-of-squares type. SPSS’s GLM procedure reports Type III by default; anova_lm’s documented default for typ is Type I. Set typ=3 and switch your categorical terms to C(x, Sum) to line the two up — changing typ alone is not enough, and will produce a third set of numbers matching neither.

My factor shows 1 degree of freedom instead of 3. What happened?

patsy read the column as continuous. This happens whenever the grouping variable is numeric — integer group codes, dose levels, timepoint numbers. Wrap it: C(group). The df column is the fastest place to catch it, before you interpret anything else.

Does typ matter for a one-way ANOVA?

No. With a single factor there is nothing to condition on, so all three types return the same F. The types only diverge with two or more factors and unbalanced cells.

Can I use anova_lm on a model I fitted with sm.OLS(y, X)?

Not reliably. anova_lm reads the model specification off the fitted result to know which design-matrix columns belong to which term, which is information only the formula interface records. Fit with statsmodels.formula.api.ols and a formula string; a hand-built design matrix will fail on the missing specification rather than give you a wrong answer.

Does statsmodels have Welch’s ANOVA?

Yes — statsmodels.stats.oneway.anova_oneway, and it is the default there (use_var='unequal'). Brown-Forsythe is use_var='bf'. The classical equal-variance test is use_var='equal'.

Does AnovaRM report Greenhouse-Geisser?

No. The docstring states that corrections for violation of sphericity are not implemented, and the output columns are only F Value, Num DF, Den DF and Pr > F. Use pingouin.rm_anova if you need epsilon and a corrected p-value.

statsmodels or pingouin?

statsmodels when the ANOVA is one output of a linear model you also want residuals, diagnostics and robust covariances from. pingouin when you want the ANOVA table a journal expects — effect sizes and sphericity corrections included — without assembling it yourself. They are not rivals: pingouin “will internally call statsmodels to calculate ANOVA with 3 or more factors, or unbalanced two-way ANOVA”.

References

• statsmodels, statsmodels.stats.anova.anova_lm — typ and robust parameters, and the Type I/II/III definitions.
• statsmodels, ANOVA user guide — the C(fcategory, Sum) worked example on the Moore dataset.
• statsmodels, statsmodels.stats.anova.AnovaRM and the anova.py source — balanced-data requirement, unimplemented between-subject effects, and the absent sphericity correction.
• statsmodels, statsmodels.stats.oneway.anova_oneway — use_var options including the Welch default.
• statsmodels, pairwise_tukeyhsd — the heteroscedasticity warning and the Games-Howell option.
• patsy, Coding categorical data — Treatment coding as the default and the C() syntax.
• Fox & Weisberg, car::Anova documentation — the contrast-coding requirement for Type III and the marginality principle behind Type II.
• SciPy, scipy.stats.f_oneway — stated assumptions and recommended alternatives.
• pingouin, anova, rm_anova and mixed_anova — ss_type, effect sizes, sphericity output and listwise deletion.
• Lakens, D. (2013). Calculating and reporting effect sizes to facilitate cumulative science. Frontiers in Psychology, 4:863.