expanding factor interactions within a formula
===

I have many formulas (of class `formula` or `Formula`) of the form `y ~ a*b`, where `a` and `b` are factors.

I need to write a function that takes such a formula and returns a formula with all of the terms in the interaction "spelled out."  Here is an example:

    fac1 <- factor('a', 'a', 'b', 'b')
    fac2 <- factor('c', 'd', 'c', 'd')
    BigFormula(formula(x ~ fac1*fac2))

where `BigFormula` returns `formula(x ~ a + b + c + d + a:c + a:d + b:c + b:d)`.

Is there a simple way to do this?

(The context: I am running many commands of the form `anova(mod1, mod2)`, where `mod2` nests in `mod1`, and where the right-hand side of both models contains terms like `fac1*fac2`.  The point of these commands is to calculate F-statistics.  The problem is that `anova` treats `fac1*fac2` as three variables when calculating degrees of freedom.  Most of the time, though, `fac1*fac2` represents more than three variables.  In the code above, for example, `fac1*fac2` represents eight variables.  As a result, `anova` is overestimating my degrees of freedom.)

