import statsmodels.api as sm from statsmodels.formula.api import ols #perform two-way ANOVA model = ols('height ~ C(water) + C(sun) + C(water):C(sun)', data=df).fit() sm.stats.anova_lm(model, typ=2) sum_sq df F PR(>F) C(water) 8.533333 1.0 16.0000 0.000527 C(sun) 24.866667 2.0 23.3125 0.000002 C(water):C(sun) 2.466667 2.0 2.3125 0.120667 Residual 12.800000 24.0 NaN NaN