To address issue #3, I have implemented now a proper formula parsing for plot2.formula:

• Internally y ~ x | z is simply converted to y ~ x + z.
• Then stats::model.frame() is called "as usual" with subset, na.action, and drop.unused.levels.
• The first variable from the resulting model frame is used as y, the second as x, and all others (if any) are collapsed to a single factor. Either via as.factor() (for a single variable) or interaction() (if more than one variable).
• Default axis and legend labels are taken from the variable names of the model frame.

Open question: Many arguments of plot2.formula() are actually not used at all in plot2.formula() but only passed on to plot2.default(). If these arguments have defaults that are identical between the two methods, the argument could also be omitted from plot2.formula(). If needed it is passed through ....