Let's start by abstracting the terminology away and look at the problem functionally in terms of f(x) = y.

x is an object containing two variables, Y_1 \& Y_2. Are they categorical or continuous? It also contains one or more variables X_i. Are the members of the set X_i all continuous, all categorical, all binary, or some combination? Some function, g measures the relationship between Y_1 \&Y_2 *or* Y_1 | Y_2 and X_i \dots X_n. As a return value y is a table of estimates and errors, etc. of main effects and interaction among factors produced by f.

Each of these roles must be mapped to the solution.

f may be `lm`

with additive (X_i + X_{i+1}) or interactive X_i * X_{i+1} terms.

g is obtained by `aov`

, which is a wrapper for f. See this helpful S/O thread

The `VGAM`

package suggestion was more for the `vignette`

. One of the difficulties in coming to `R`

from SAS or another platform is mapping from the familiar to the unfamiliar. Approaching the task as function composition of objects to objects helps clarify the syntactical requirements in the new language.