R: new formula & vegan breaks

[jarioksa]

I just got the following email message from Martin Maechler:

just a quick note before I embark on longish trip to Calgary Canada:
Following up the thread on the R-devel mailing list,
I've recently changed formula(<model.frame>) behavior such that it
returns the formula from the "terms" attribute of the model frame,
rather than a newly constructed formula from the model.frame's
variables (as it does for non-model.frame data frames).

You may have noticed on the CRAN checks for vegan that this breaks at
least one of vegan's examples on "R-devel" (R's development version).

I may take time on the long plane flight to see how to fix it... or
may not: vegan seems to be doing pretty "strange" contortions in
evaluating formulas,
as far as I could tell, quite a bit non-standard and different than
say lm(), glm() etc.

This may mean that we have to do something with the formula. Most likely this concerns cca.formula and its ilk. Other functions have more usual handling of formula, although we do have matrix-like right-hand-side in most functions.

[gavinsimpson]

The proximal issue of the error is in ordiParseFormula() on L76. There we call

mf <- model.frame(formula(mf), mf, xlev = xlev,
                  na.action = na.action, drop.unused.levels = TRUE)

and in R-patched formula(mf) returns

> formula(mf)
Management ~ `poly(A1, 2)` + spno

whereas in R-devel it returns

> formula(mf)
~Management + poly(A1, 2) + spno

The current behaviour seems wrong — why is Management now the response? — but the direct issue is that currently we expect poly(A1, 2) to be quoted as we want it to find the variable (a matrix) of that same name in mf:

> str(mf)
'data.frame':   13 obs. of  3 variables:
 $ Management : Factor w/ 4 levels "BF","HF","NM",..: 1 4 4 NA 2 2 1 4 4 3 ...
 $ poly(A1, 2): num [1:13, 1:2] -0.1421 -0.0579 -0.0684 0.1526 -0.0579 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr  "1" "2"
....

instead of re-evaluating the formula — we don't store A1 or at least it's not available by the time we evaluate the model frame on L76. If it was R-devel would work.

[jarioksa]

There is a plenty of cruft (and nothing's plenty for me) in ordiParseFormula. Some of this is just to get around quirks in R formula and model.frame behaviour when we need to pass the structures to various functions. For instance, capscale was calling rda.default in previous versions of vegan so that the data was handled in different environment (depth) in rda.default when it was called from rda directly or via capscale. Things have been simplified a lot and the interface unified in the current vegan allowing for simplification (but the old code still worked and was not changed).

About the formula(mf) issue: this is a change in R behaviour, and I am not sure if it is important for our case. The current R just assumes that model frame has a response, and takes the first variable to be that response. We removed the response (community data), but current R (release) does not know about this possibility. The original formula without the lhs can be extracted with

## R 3.5.x
formula(terms(mf))

but stats:::formula.data.frame does not query this, but it reconstructs the formula from names(mf), always taking the first name as the lhs, and quoting the names like poly(A1, 2). So I guess that R-devel has switched to different way of extracting the formula.

[jarioksa]

This really seems to be an issue with changes in stats:::formula.data.frame behaviour. Earlier it reconstructed the formula from model.frame names and back-quoted names which contained functions, but now it just fetches the saved formula. With current scoping in vegan, the new behaviour fails. This is pretty easy to fix by adding still one kludge to vegan:::ordiParseFormula, and @mmaechler even provides compatibility function stats::DF2formula to keep the old behaviour. However, I now suggest a fix that explicitly takes care of having back-quotes when we need them, and we do not need to query R version.

The changes are in branch https://github.com/vegandevs/vegan/tree/issue-%23299-kludge

This is yet another kludge in a kludgy function and the real solution is to re-write (or scrap) vegan:::ordiParseFormula. This is not something that you can do while having your afternoon coffee break, but needs time and will break vegan for a long time before getting everything to work again. Someday, someone, somehow...

[gavinsimpson]

This looks good @jarioksa I mean ideally we'd rewrite it but that is a longer-term project and if this maintains current behaviour then that's good.

[jarioksa]

Closed with commit 0e7b8be. This is a kluge, and in long term, whole ordiParseFormula should be re-designed and re-written.