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.