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