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