2
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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])
$\endgroup$
  • $\begingroup$ Thanks - turns out that I used the wrong contrasts. I guess they weren't orthogonal so I didn't get the correct SSs. $\endgroup$ – cbrnr Jul 25 '16 at 12:22
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ This makes sense. However, how do I calculate the interaction SS then? In the two-way case it's obviously just the full model minus the two main effects - but how does this generalize to the N-way case? If I subtract all N main effects from the full model, I end up with all possible interactions. $\endgroup$ – cbrnr Jul 18 '16 at 6:58
  • $\begingroup$ For a balanced design you can do it all by hand but it is tiresome. For an unbalanced design it is worse. Using a computer is really the way to go here. $\endgroup$ – mdewey Jul 18 '16 at 18:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.