expanding factor interactions within a formula

Look at the help for formula there may be existing things that will work for you.

For example the formula y ~ (a + b + c + d)^2 will give you all main effects and all 2 way interactions and the formula y ~ (a + b) * (c + d) gives the expansion that you show above. You can also subtract terms so y ~ a*b*c - a:b:c will not include the 3 way interaction.

How about the following solution. I use a more extreme example of a complex interaction.

f = formula(y ~ a * b * c * d * e)

To spell out the interaction terms, we extract the terms from the value returned by terms.formula():

terms = attr(terms.formula(f), "term.labels")

which yields:

> terms
 [1] "a"         "b"         "c"         "d"         "e"         "a:b"       "a:c"      
 [8] "b:c"       "a:d"       "b:d"       "c:d"       "a:e"       "b:e"       "c:e"      
[15] "d:e"       "a:b:c"     "a:b:d"     "a:c:d"     "b:c:d"     "a:b:e"     "a:c:e"    
[22] "b:c:e"     "a:d:e"     "b:d:e"     "c:d:e"     "a:b:c:d"   "a:b:c:e"   "a:b:d:e"  
[29] "a:c:d:e"   "b:c:d:e"   "a:b:c:d:e"

And then we can convert it back to a formula:

f = as.formula(sprintf("y ~ %s", paste(terms, collapse="+")))

> f
y ~ a + b + c + d + e + a:b + a:c + b:c + a:d + b:d + c:d + a:e + 
    b:e + c:e + d:e + a:b:c + a:b:d + a:c:d + b:c:d + a:b:e + 
    a:c:e + b:c:e + a:d:e + b:d:e + c:d:e + a:b:c:d + a:b:c:e + 
    a:b:d:e + a:c:d:e + b:c:d:e + a:b:c:d:e