Calculate interaction sum of squares by hand I'm trying to manually calculate the sums of squares (SSs) in a two-way ANOVA in Python. However, I cannot seem to get the correct result for the interaction SS. I'm using the ToothGrowth dataset as an example, which has two factors $A$ (supp, two levels) and $B$ (dose, three levels).
Dummy coding $A$ requires one dummy variable, whereas $B$ requires two dummy variables. Coding the interaction involves multiplying the dummies for $A$ and $B$, thus yielding two new dummy variables.
The full model is $SS_F = SS_A + SS_B + SS_{AB}$. It consists of all five dummy variables. When I compute the SSs for the full model and for the three reduced models, I get the following result:

              ssm          ssr          sst
full  2740.103333   712.106000  3452.209333
supp   205.350000  3246.859333  3452.209333
dose  2426.434333  1025.775000  3452.209333
ab     648.336583  2803.872750  3452.209333

Obviously, the interaction SS is wrong, because it should be equal to $SS_F - SS_A - SS_B = 108.319$. In my calculations, I use the two dummy predictors for the interaction (and an intercept) to get the interaction SS. What's wrong with this - why is the full model involving the same interaction variables correct but the interaction only is wrong?
 A: An interaction means that the effect of dose depends on what level of supp is applied with it. If you just include the two dummies which correspond to the interaction in the full model you are not estimating that. It might be possible to attach some meaning to your model although I am not sure what that might be but it is not interaction in the usual sense of that word.
A: You must have made a mistake in calculating the interaction sum of squares. I am by no means an expert on this topic, but the following Python code gives the expected results:
import numpy as np
import pandas as pd


def factor_sos(A, coef, idx):
    factor_contribution = np.dot(A[:, idx], coef[idx])
    sos_explained = np.sum(factor_contribution**2)
    return sos_explained    


data = pd.read_csv('~/Downloads/ToothGrowth.csv')
dummy = pd.get_dummies(data, columns=['dose', 'supp'])

# Contrasts
supp = dummy[['supp_OJ']].values - dummy[['supp_VC']].values
dose = dummy[['dose_1.0', 'dose_2.0']].values - dummy[['dose_0.5']].values
interaction = supp * dose

A = np.hstack([supp, dose, interaction])
b = dummy['len'].values
coef, res, _, _ = np.linalg.lstsq(A, b)

sos_full = factor_sos(A, coef, [0, 1, 2, 3, 4])
sos_supp = factor_sos(A, coef, 0)
sos_dose = factor_sos(A, coef, [1, 2])
sos_interaction = factor_sos(A, coef, [3, 4])

print('Sum of squares:')
print('full      ', sos_full)
print('supp      ', sos_supp)
print('dose      ', sos_dose)
print('supp*dose ', sos_interaction)

# Sum of squares:
# full       2740.10333333
# supp       205.35
# dose       2426.43433333
# supp*dose  108.319

First, this code fits a linear regression model that includes all dummy variables. Then, the function factor_sos is used to calculate how much each factor and interaction of factors contributes to the variability (sum of squares) of the model's output.
If you are only interested in a specific factor or interaction, there is no need to fit the full model. This works too:
A = np.hstack([interaction])
coef, res, _, _ = np.linalg.lstsq(A, b)
sos_interaction = factor_sos(A, coef, [0, 1])

